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ABSTRACT 

Current theories predict relativistic hadronic particle populations in clusters of galaxies in addition to the already 
observed relativistic leptons. In these scenarios hadronic interactions give rise to neutral pions which decay into 
y rays that are potentially observable with the Large Area Telescope (LAT) on board the Fermi space telescope. 
We present a joint likelihood analysis searching for spatially extended ]/-ray emission at the locations of 50 galaxy 
clusters in four years of Fermi-LAI data under the assumption of the universal cosmic-ray (CR) model proposed 
by Pinzke & Pfrommer. We find an excess at a significance of 2.7a, which upon closer inspection, however, 
is correlated to individual excess emission toward three galaxy clusters: A400, A1367, and A3 112. We discuss 
these cases in detail and conservatively attribute the emission to unmodeled background systems (for example, 
radio galaxies within the clusters). Through the combined analysis of 50 clusters, we exclude hadronic injection 
efficiencies in simple hadronic models above 21% and establish limits on the CR to thermal pressure ratio within 
the virial radius, R 200 , to be below 1.25%-1.4% depending on the morphological classification. In addition, we 
derive new limits on the y-ray flux from individual clusters in our sample. 

Key words: galaxies: clusters: intracluster medium - gamma rays: galaxies: clusters 
Online-only material: color figures 


1. INTRODUCTION 

The quest for the first detection of high-energy y rays from 
galaxy clusters is still ongoing. While there have been y-ray 
detections from radio galaxies in clusters such as NGC 1275 
(Strong & Bignami 1983; Abdo et al. 2009; Aleksic et al. 2012b) 
and IC 310 (Aleksic et al. 2010a; Neronov et al. 2010) in the 
Perseus cluster, as well as M87 in the Virgo cluster (Sreekumar 
et al. 1994; Abdo et al. 2009; Aharonian et al. 2003), no cluster- 
wide y-ray emission has been detected so far. Previous reports 
of space-based cluster observations in the GeV band include 
Reimer et al. 2003, Ackermann et al. 2010b, Zimmer et al. 
2011, Ando & Nagai 2012, and Han et al. 2012. Ground-based 
observations in the energy band >100 GeV were reported for 
Perseus and A2029 by Perkins et al. (2006), for A496 and A85 
by Aharonian et al. (2009b), for Coma by Aharonian et al. 
(2009a); Arlen et al. (2012), for A3667 and A4038 by Kiuchi 
et al. (2009), for Perseus by Aleksic et al. (2010a, 2012a), 
and for Fornax by Abramowski et al. (2012). Nevertheless, 
current space-based y-ray detectors such as the Large Area 
Telescope (LAT) on board the Fermi satellite (Atwood et al. 
2009) may be able to detect y rays from galaxy clusters during its 
lifetime. 

The discovery and characterization of cosmic-ray-induced 
(CR-induced) y rays from clusters could not only serve as 
a crucial discriminator between different models for the ob- 
served cluster- wide radio emission (Pfrommer & EnBlin 2004a; 
Brunetti et al. 2012) but could also be a signpost of the physical 
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heating processes underlying feedback by active galactic nuclei 
(AGNs) that has been proposed to solve the “cluster cooling flow 
problem” (Loewenstein et al. 1991; Guo & Oh 2008; EnBlin et al. 
201 1 ; Fujita & Ohira 2012; Wiener et al. 2013; Pfrommer 2013). 
Moreover, any observed y-ray emission from clusters necessar- 
ily requires a detailed understanding of the astrophysical y-ray 
contribution before putting forward constraints on more exotic 
scenarios such as annihilation or decay signals from dark matter 
(DM; Ackermann et al. 2010a; Dugger et al. 2010; Pinzke et al. 
2011; Huang et al. 2012). 

The intracluster medium (ICM) consists of a hot (1-10 keV) 
plasma, which has primarily been heated through collisionless 
shocks that form as a result of the hierarchical build-up of galaxy 
clusters (e.g., Ryu et al. 2003; Pfrommer et al. 2006; Skill- 
man et al. 2008; Vazza et al. 2009). These structure formation 
shocks and the turbulent motions of the gas in combination 
with intracluster magnetic fields provide the necessary condi- 
tions for efficient particle acceleration (e.g., Colafrancesco & 
Blasi 1998; Ryu et al. 2003). Recent radio observations prove 
the existence of relativistic electrons and magnetic fields in 
clusters. Some cool core (CC) clusters host radio mini halos 
in their centers (EnBlin et al. 2011). Additionally, a subsample 
of merging non-cool core (NCC) clusters show radio relics at 
their periphery (Kempner et al. 2004) and/or giant radio halos 
that often extend out to Mpc scales. While giant radio halos 
have been observed in more than 50 clusters (see, e.g., Ferrari 
et al. 2008, for a review), their precise origin is still not under- 
stood. There are two competing theories to explain these radio 
halos. 

In hadronic models, CR ions and protons {p) are acceler- 
ated in structure formation shocks, jets of radio galaxies, and 
supernovae-driven galactic winds (see, e.g., Volk et al. 1996; 
EnBlin et al. 1997; Berezinsky et al. 1997; Pfrommer & EnBlin 
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2004a), and significant populations of CR protons can accumu- 
late due to their long cooling time in the ICM (Volk et al. 1996). 
Inelastic collisions of CR ions with thermal protons of the ICM 
produce both neutral and charged pions, which decay almost 
instantly into y rays and electrons/positrons, respectively. This 
process could in principle account for the radio-emitting leptons, 
while requiring only a modest CR-to-thermal pressure ratio of 
(at most) a few percent (Dennison 1980; Vestrand 1982; Blasi 
& Colafrancesco 1999; Dolag & EnBlin 2000; Miniati et al. 
2001a, 2001b; Miniati 2003; Pfrommer & EnBlin 2003, 2004a, 
2004b; Blasi et al. 2007; Pfrommer et al. 2008; Pfrommer 2008; 
Kushnir et al. 2009; Donnert et al. 2010a, 2010b; Keshet & Loeb 
2010; Keshet 2010; EnBlin et al. 2011). The non-detection of 
y-ray emission from individual radio halo clusters places strong 
limits on intracluster magnetic fields within the hadronic model 
(Pfrommer & EnBlin 2004b; Jeltema & Profumo 2011; Arlen 
et al. 2012; Brunetti et al. 2012; Aleksic et al. 2012a). The sec- 
ondary leptons also experience Compton upscattering with the 
background radiation fields to y-ray energies, but this expected 
emission is always subdominant compared to the y rays pro- 
duced by decaying neutral pions (Miniati 2003; Pfrommer et al. 
2008; Pinzke & Pfrommer 2010). 

The re-acceleration models assume the existence of a long- 
lived pool of mildly relativistic electrons, that were accelerated 
in the past by structure formation shocks, galactic winds, and 
AGNs, or coincide with secondary electrons that are injected 
in the aforementioned hadronic CR p-p interactions. Those CR 
electrons scatter with plasma waves that are excited by ICM 
turbulence, e.g., after a cluster merger. These particle-wave 
interactions may accelerate the particles through the second 
order Eermi process to sufficiently high energies to explain the 
observed radio emission (Schlickeiser et al. 1987; Brunetti et al. 
2001; Petrosian 2001; Brunetti et al. 2004; Brunetti & Blasi 
2005; Brunetti & Lazarian 2007, 2011; Brunetti et al. 2009; 
Donnert et al. 2013). 

Assuming that the same physical processes that produce y 
rays are present in each galaxy cluster, independent of mass, 
age, and other characteristics, we employ a joint likelihood 
analysis to search for these y rays. The resulting universal 
scaling factor Ay can be used to derive limits on the hadronic 
acceleration efficiency at structure formation shocks and the 
volume-averaged CR-to-thermal pressure (Xcr). While the 
joint likelihood method can be applied to study any emission 
governed by a universal physical process, in this paper, we focus 
on the search for y rays from CR-induced pion decay and defer 
more exotic scenarios such as y rays from DM annihilation and 
decay to future studies. 

The outline of this paper is as follows. We discuss our cluster 
selection in Section 2 and address the Fermi-hAI observations 
and data analysis in Section 3. Our cluster emission models are 
described in Section 4. We present and discuss our results in 
Section 5. Einally, we present a survey of possible systematics 
in Section 6, and we conclude in Section 7. Throughout the 
paper, we assume a cosmology with Qm = 0.3, Qa = 0.7, and 
h = 0.7. 

2. CLUSTER SELECTION 

Assuming a correlation between y-ray and X-ray luminosity 
in galaxy clusters (Colafrancesco & Blasi 1998; Pfrommer & 
EnBlin 2004a), we use the extended Highest X-ray Elux Galaxy 
Cluster Sample (HIELUGCS; Reiprich & Bohringer 2002; Chen 
et al. 2007), containing the 106 nearby brightest X-ray clusters 


from the ROSAT all- sky survey, as a suitable list of sources when 
constructing our analysis sample. 

Given that our statistical approach assumes independent 
sources, we remove clusters where the angular separation 
between cluster centers is less than any of the respective virial 
radii enlarged by 1°.^^ This cut is motivated by Monte Carlo 
(MC) studies (see Appendix A. 1 for details), showing that for an 
energy threshold of 500 MeV, the bias on the likelihood ratio test 
statistic due to overlaps is minimal under the condition above. 
However, in the case that the expected y-ray flux from a single 
HIELUGCS cluster within such an ensemble of overlapping 
clusters is responsible for more than 90% of the total expected 
emission, as calculated using the approach in Pinzke et al. 
(2011), we neglect the other clusters in the ensemble and 
attribute all photons to the cluster with the largest expected flux. 

The aforementioned virial radius, R 200 , of the cluster is 
the radius containing the virial mass, M 200 , which in turn 
is derived from the M 500 mass reported by Chen et al. 
(2007).^^ We solve for M 200 using M 200 = ^200 x 200/500 x 
[c 2 oo(^ 20 o)/^^ 50 o(^ 20 o)]^ (Voit 2005), where the halo-mass- 
dependent concentration parameter c is derived from a power- 
law fit to a sample of observed galaxy cluster concentration and 
masses (Comerford & Natarajan 2007). 

The extended HIELUGCS catalog contains clusters up to 
redshifts of 0.18, with the majority located at z < 0.1. Eor 
simplicity, we do not apply a redshift correction to the spectrum, 
and we decided to exclude all clusters above z = 0.1, as 
their inclusion without applying a correction would amount to 
incorrect modeling of the spectrum. These clusters contribute 
less than 1 % to the total expected y -ray flux including all clusters 
in the HIELUGCS catalog, as derived in Pinzke et al. (2011), 
making this a suitable approach. In addition, we exclude clusters 
that lie in a box defined by | Z? | ^50° and | / 1 ^ 20° , as this region 
contains emission from the Galactic plane and the Fermi bubbles 
(Su et al. 2010), the models for which have relatively large 
systematic uncertainties. The Virgo cluster, which has a virial 
radius of ~8° on the sky (Eouque et al. 2001) and clusters that 
fall into this region, such as M49, is excluded from the analysis. 

The Galactic plane presents a substantial challenge due 
to a large number of potentially unresolved point sources 
and uncertainties in modeling the diffuse foregrounds. We 
conservatively define the plane region to be \b\ ^ 20° and 
remove clusters from our sample that lie within this region, such 
as the Centaurus cluster. There are nine Abell clusters present 
in the HIELUGCS catalog that are located within a radius of 7° 
from the center of the Centaurus cluster which overlap with one 
another. In order to avoid highly crowded analysis regions, we 
exclude the entire region. 

The remaining clusters are considered separately according to 
their classification as either CC, NCC, or unclassified. Eor this 
purpose, we use the classification presented in Hudson et al. 
(2010) or follow recommendations by Cavagnolo et al. (2009). 
Among other quantities, we classify clusters that have central 
entropy values Kq ^ 30keV cm^ as CC and otherwise as NCC. 
When there is no supporting X-ray data available, we leave these 
clusters unclassified. 

The 50 galaxy clusters we consider in this analysis are 
listed in Table 1 along with characteristic cluster quantities 


All cluster positions, which were taken from the NASA Extragalactic 
Database (http://ned.ipac.caltech.edu/), are based on observations in the optical 
waveband. 

We define the virial radius of a cluster as the radius at which the mean 
interior density equals 200 times the critical density of the universe. 
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Table 1 

Cluster Sample Considered for the Analysis 


Name 

R.A. 

(°) 

Deck 

(°) 

z 

R 200 

(Mpc) 

M 200 

(xlO‘5) 

Ext. 

(°) 

FCRp(£ > 500 MeV) 
(10-'Vcm-2 s“') 

(^CR>«hl 

(xl0-‘) 

(^CR>fi200 

(xl0-‘) 

Morph. 

2A0335 

54.65 

9.97 

0.04 

1.31 

0.26 

1.06 

2.84 

0.32 

0.14 

CC 

A0085 

10.41 

-9.34 

0.06 

1.89 

0.78 

1.00 

2.93 

0.24 

0.10 

cc 

A0119 

14.09 

-1.26 

0.04 

1.96 

0.86 

1.27 

1.42 

0.24 

0.14 

NCC 

A0133 

15.66 

-21.96 

0.06 

1.52 

0.41 

0.78 

0.72 

0.28 

0.13 

CC 

A0262 

28.21 

36.15 

0.02 

0.91 

0.09 

1.54 

1.21 

0.41 

0.30 

CC 

A0400 

44.41 

6.03 

0.02 

1.02 

0.12 

1.17 

0.44 

0.38 

0.28 

NCC 

A0478 

63.34 

10.48 

0.09 

1.95 

0.86 

0.67 

2.45 

0.24 

0.09 

CC 

A0496 

68.41 

-13.25 

0.03 

1.58 

0.46 

1.36 

2.83 

0.28 

0.12 

CC 

A0548e 

87.16 

-25.47 

0.04 

1.06 

0.14 

0.76 

0.25 

0.37 

0.26 

7 

A0576 

110.35 

55.74 

0.04 

1.56 

0.44 

1.14 

0.63 

0.28 

0.16 

NCC^ 

A0754 

137.21 

-9.64 

0.05 

2.27 

1.35 

1.22 

1.69 

0.21 

0.08 

NCC 

A1060 (Hydra) 

159.21 

-27.53 

0.01 

1.26 

0.23 

2.77 

2.28 

0.33 

0.18 

CC 

A1367 (Leo) 

176.12 

19.84 

0.02 

1.83 

0.71 

2.33 

1.92 

0.25 

0.13 

NCC 

A1644 

194.31 

-17.35 

0.05 

1.83 

0.70 

1.11 

1.30 

0.25 

0.14 

CC 

A1795 

207.25 

26.59 

0.06 

2.02 

0.95 

0.95 

3.01 

0.23 

0.11 

CC 

A2065 

230.68 

27.72 

0.07 

2.11 

1.09 

0.86 

0.92 

0.21 

0.08 

CC 

A2142 

239.57 

27.22 

0.09 

2.30 

1.40 

0.77 

3.45 

0.29 

0.14 

CC 

A2199 

247.16 

39.55 

0.03 

1.52 

0.40 

1.42 

3.00 

0.27 

0.13 

CC 

A2244 

255.68 

34.05 

0.10 

1.66 

0.52 

0.52 

0.76 

0.25 

0.14 

CC 

A2255 

258.13 

64.09 

0.08 

1.87 

0.76 

0.69 

0.85 

0.22 

0.11 

NCC^ 

A2256 

255.93 

78.72 

0.06 

2.17 

1.18 

1.09 

2.55 

0.31 

0.17 

NCC^ 

A2589 

351.00 

16.82 

0.04 

1.38 

0.30 

0.95 

0.68 

0.30 

0.13 

CC 

A2597 

351.33 

-12.11 

0.09 

1.45 

0.35 

0.51 

0.76 

0.28 

0.17 

CC 

A2634 

354.58 

27.03 

0.03 

1.55 

0.43 

1.39 

0.62 

0.26 

0.13 

CC 

A2657 

356.21 

9.14 

0.04 

1.71 

0.58 

1.21 

0.83 

0.28 

0.15 

CC 

A2734 

2.83 

-28.87 

0.06 

1.58 

0.46 

0.74 

0.44 

0.26 

0.12 

NCC 

A2877 

17.46 

-45.90 

0.03 

1.78 

0.65 

2.03 

0.51 

0.28 

0.12 

NCC^ 

A3112 

49.47 

-44.24 

0.08 

1.53 

0.41 

0.60 

1.11 

0.27 

0.14 

CC 

A3158 

55.67 

-53.63 

0.06 

1.68 

0.55 

0.82 

1.20 

0.20 

0.09 

NCC 

A3266 

67.80 

-61.41 

0.06 

2.54 

1.89 

1.26 

3.13 

0.26 

0.15 

CC 

A3376 

90.18 

-40.05 

0.05 

1.78 

0.65 

1.12 

0.69 

0.27 

0.17 

NCC 

A3822 

328.53 

-57.85 

0.08 

1.57 

0.45 

0.61 

0.45 

0.28 

0.17 

NCC^ 

A3827 

330.45 

-59.95 

0.10 

2.36 

1.52 

0.73 

0.98 

0.21 

0.09 

NCC^ 

A3921 

342.41 

-64.39 

0.09 

1.77 

0.63 

0.58 

0.47 

0.26 

0.13 

NCC^ 

A4038 

356.90 

-28.13 

0.03 

1.28 

0.24 

1.20 

1.34 

0.32 

0.17 

CC 

A4059 

359.17 

-34.67 

0.05 

1.54 

0.42 

0.93 

0.97 

0.28 

0.13 

CC 

COMA 

194.95 

27.98 

0.02 

2.02 

0.96 

2.46 

11.40 

0.23 

0.12 

NCC 

EXO0422 

62.93 

-29.81 

0.04 

1.30 

0.25 

0.93 

0.72 

0.32 

0.17 

CC 

EORNAX 

54.63 

-35.46 

0.01 

1.01 

0.12 

5.98 

0.84 

0.38 

0.25 

CC 

HCG94 

349.32 

18.72 

0.04 

1.22 

0.21 

0.83 

0.38 

0.33 

0.21 

7 

HYDRA-A 

139.52 

-12.10 

0.06 

1.50 

0.38 

0.79 

1.67 

0.29 

0.12 

CC 

IIIZw54 

55.32 

15.40 

0.03 

1.45 

0.35 

1.41 

0.45 

0.30 

0.16 

cc 

IIZwl08 

318.48 

2.57 

0.05 

1.47 

0.36 

0.86 

0.40 

0.29 

0.19 

7 

NGC 1550 

64.91 

2.41 

0.01 

0.81 

0.06 

1.80 

0.67 

0.45 

0.33 

CC 

NGC 5044 

198.85 

-16.39 

0.01 

0.73 

0.04 

2.14 

0.61 

0.50 

0.42 

cc 

RXJ2344 

356.07 

-4.37 

0.08 

1.95 

0.86 

0.74 

0.57 

0.24 

0.11 

cc 

S405 

57.89 

-82.22 

0.06 

1.56 

0.44 

0.74 

0.40 

0.28 

0.18 

cc^ 

S540 

85.03 

-40.84 

0.04 

1.27 

0.24 

1.00 

0.35 

0.32 

0.18 

NCC 

UGC03957 

115.24 

55.43 

0.03 

1.39 

0.31 

1.15 

0.50 

0.30 

0.15 

7 

ZwC11742 

266.06 

32.98 

0.08 

2.04 

0.98 

0.80 

1.17 

0.23 

0.10 

CC 


Notes. The sample we used to carry out our analysis. The locations are taken from the NED database. The columns from left to right read: cluster name 
according to Reiprich & Bohringer (2002); right ascension (J2000); declination (J2000); redshift, R 200 ', M 200 in units of Mq/H^q; extension (2 x ^200, 
where ^200 refers to the angular virial radius as described in footnote 62); expected y-ray flux above 500 MeV; CR-to-thermal pressure ratio within 
R 200 and the half-light radius, Rul (see Section 4.1 and footnote 66 for details); and morphological classiflcation. We denote unclassifled clusters 
with a question mark. We note that columns 8-10 are based on the predictions in Pinzke et al. (2011) and assume Ay = 1. Redshifts are from Chen 
et al. (2007) and references therein. M200 and R 200 are calculated from R500 and M500 in the aforementioned reference. Unless specifled otherwise, all 
classiflcations are from Hudson et al. (2010). 

^ Values and classiflcations are from Sivanandam et al. (2009). 

^ Values and classiflcations are taken from the ACCEPT database (Cavagnolo et al. 2009). 
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Figure 1. Summary of cluster quantities in our analysis. In the upper panel, we show redshift vs. M 200 and in the lower panel the predicted y-ray flux above 500 MeV 
taken from Pinzke et al. (2011), > 500 MeV) vs. its extension (= 2 x ^ 200 ) in degrees. Most of the selected clusters in the sample have extensions of ~1°. 



Figure 2. Hammer- Aitoff projection of the sky as seen by the LAT after four years of exposure. Shown is a counts map generated for CLEAN class events in the energy 
range from 500 MeV to 200 GeV using the same set of quality cuts as described in Section 3. The dashed circles correspond to the analysis regions considered for 
this analysis. The solid circles represent the clusters used in this analysis and their extension as characterized by their virial radii. In red we show CC, in green NCC, 
and in blue unclassifled clusters. We shade the region which is described by the Fermi bubbles in Su et al. (2010) while schematically overlaying our geometric cuts 
for masking the Galactic plane (\b\ ^ 20°) and the Galactic bubbles (\b\ ^ 50° and |/| ^ 20°). The latter has been designed to mask out the majority of the emission 
that can be attributed to the lobes. When comparing our mask with the geometric description by Su et al. (2010), we found that two clusters are contained inside the 
bubbles. We checked that the region containing these two clusters is modeled well and hence keep the clusters in our analysis. 

(A color version of this flgure is available in the online journal.) 


(see Figure 1) and their classification. We show their location 
and radial extension on the sky in Figure 2.^^ 


The extension we use is twice the angle subtended by the angular virial 
radius, O 200 , which is given by O 200 = arctan(/? 20 o/A>a) x 180 ° /tt, where Da is 
the angular diameter distance from the Earth to the center of the cluster. 


3. Fermi-LAT OBSERVATIONS AND DATA ANALYSIS 

Fermi-LAI is the main instrument on board the Fermi 
Gamma-ray Space Telescope. Since its launch in 2008, the LAT 
has surveyed the y-ray sky in the energy range from 20 MeV 
to >300 GeV with unprecedented sensitivity. For more details 
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about the LAT, the reader is referred to Atwood et al. (2009) and 
to Ackermann et al. (2012a) for the on-orbit LAT performance. 
We carry out a binned likelihood analysis of 48 months of 
Fermi-LKV data (2008 August 4-2012 August 4), which has 
been reprocessed to account for the time-dependent calorimeter 
response (Bregeon et al. 2013), and we refer to this data as 
P7REP. We select events corresponding to the CLEAN class, 
which consists of the events that have the highest probability of 
being y rays, and we use the P7REP_V15 instrument response 
functions (IRFs) provided in the software package Fermi 
Science Tools v9r32p5.^^ We apply standard quality cuts 
to our data using gtmktime and require DATA_QUAL==1 && 
LAT_C0NFIG==1, which refers to the configuration during 
nominal science operation. We require the magnitude of the 
rocking angle of the LAT to be ^52° and reject events above a 
zenith angle of 100° to greatly reduce contamination by Earth 
limb emission. We use gtbin to bin the data in 0? 1 spatial bins.^"^ 
The spectra are binned in 1 8 logarithmically spaced bins from 
500 MeV to 200 GeV. Above 200 GeV, the number of expected 
events given the models and the number of detected events are 
both sufficiently low that the models are not well constrained, 
so we omit this energy range from our analysis. 

3.1. Joint Likelihood 

The joint likelihood is a source stacking technique, which 
has been previously applied with LAT data in the search for 
DM (Ackermann et al. 2011) and to study the extragalactic 
background light (Ackermann et al. 2012b). In brief, if the goal 
is to constrain or estimate a single or a set of parameters common 
to a source class, then backgrounds and individual properties of 
each source can be modeled individually and treated as nuisance 
parameters in source- specific likelihoods. The source-specific 
likelihoods can then be multiplied to yield a joint likelihood 
function that is used for inference on the common parameter of 
interest. The joint likelihood function for our case can be written 
as 

C{Ay\V) = Y\Ci{Ay,h\Vi), (1) 

i 

where Ay is a. (dimensionless) universal scale factor which 
serves as the parameter of interest and V refers to the photon data 
for all regions-of-interest (ROIs). The physical interpretation of 
the universal scale factor in terms of CR-induced y-ray emission 
is discussed further in Section 4.1. The bi parameters correspond 
to the parameters describing the background components in 
the individual ROIs and are treated as nuisance parameters in 
the likelihood evaluation. The hiS include the normalizations 
of the isotropic and Galactic diffuse components. We denote 
the photon data for each ROI as Vi . The index i runs over the 
ROIs in the sample. Having constructed the likelihood function, 
we use the profile likelihood method (e.g., Rolke et al. 2005) 
to obtain the best-fit values and confidence intervals for the 
parameter of interest. Ay . We constrain our parameter of interest 
and nuisance parameters to be positive. The joint likelihood 
function is implemented in the Fermi Science Tools using 
the Compos it e2 package, and profiling over the likelihood 
function is achieved by means of MINOS, which is part of the 


MINUIT package (James & Roos 1975). The 90% confidence 
interval, [A^, is defined as the values of Ay where the 
log-likelihood has changed by 2.71/2 with respect to its value 
at the maximum, — 2Alog C = 2.1 \ (Bartlett 1953; Rolke et al. 
2005). These intervals can be reinterpreted as upper limits at the 
95% confidence level (C.L.), if the parameter is unconstrained 
in the fit, which we do if the lower limit ^0. 

For quantifying the significance of a potential excess, we 
employ the common likelihood ratio approach (Neyman & 
Pearson 1928): 


TS = -2 log 


/ C{Ay = 0 , 1) 

\ C{Ay^b) 


( 2 ) 


where C{Ay = 0, b) is the null hypothesis, i.e., it represents 

the likelihood evaluated at the best-fit b under the background- 
only hypothesis (Section 5.1) and C{Ay,b) is the likelihood 
evaluated at the best-fit value of Ay, b, when including our 
candidate y-ray (cluster) source. 

As we constrain the signal fit parameter Ay ^ 0, the null 
distribution of the test statistic, TS, is given by (1/2)5 -i- (l/2)x^ 
for one degree of freedom (Chemoff 1952), which we veri- 
fied by MC simulations. For details, the reader is referred to 
Appendix A. While these simulations agree well with the ex- 
pectations from Chernoff’s theorem, studies based on random 
ROIs that encapsulate systematic effects in the LAT data (e.g., 
imperfect diffuse modeling, unresolved background sources, 
and percent-level inconsistencies in the IRFs) indicate that 
the significance estimated by simulations is probably some- 
what too high when compared to the asymptotic expectations 
(Ackermann et al. 2014). 

Recalling the discussion in Section 2, sources and ROIs have 
to be selected such that the overlap is minimal since the joint 
likelihood function Equation (1) does not account for correlation 
terms between the different individual likelihood functions. For 
technical reasons, we cannot define a single ROI containing 
all 50 clusters, as this leads to an overflow in the number of 
free parameters that MINOS is not able to handle. We therefore 
construct non-overlapping ROIs, each containing one or more 
(non-overlapping) cluster sources. 

3.2. Construction of Regions of Interest 

In order to avoid overlaps between ROIs, we perform an 
iterative procedure in which we treat each cluster in our sample 
with its extension as listed in Table 1 as a seed source and 
construct a circular region around it. We require these regions 
to be at least 8° in radius in order to obtain a good fit to 
the background model. If we cannot accommodate clusters in 
separate ROIs, we mitigate any remaining overlap by enlarging 
the ROIs such that more than one cluster may be contained in 
one analysis region. 

Using this method, we are able to construct 26 independent 
ROIs containing the clusters in our sample which are listed in 
Table 2. We show the locations and extensions of our selected 
clusters and the ROIs containing them in Figure 2. 


Both the LAT data as well as the appropriate analysis tools are made 
available to the public by the Fermi Science Support Center 
http://fermi.gsfc.nasa.gov/ssc/data/. 

It should be noted that the binned likelihood analysis using the Fermi 
Science Tools uses square shaped analysis regions. For simplicity, we refer 
to our ROIs by the central coordinates and radii of the ROI circumcircles. 


4. SIGNAL AND BACKGROUND MODEL 

4.1. Signal Model: Gamma-Ray Emission from Cosmic Rays 

In this paper, we focus on the CR-induced y-ray signal 
and defer an analysis of y rays originating from DM decay 
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Table 2 

Regions of Interest Considered in This Analysis 


Region 

R.A. 

(°) 

Decl. 

(°) 

Radius 

o 

Clusters 

cc 

NCC 

1 

12 

-5 

8 

2 

1 

1 

2 

15.66 

-21.95 

8 

1 

1 


3 

354 

-9 

8 

2 

2 


4 

172.07 

21.02 

8 

1 


1 

5 

30.3 

35.75 

8 

1 

1 


6 

257.77 

73.43 

16 

2 


2 

7 

235 

27.5 

8 

2 

2 


8 

259.5 

40 

16 

3 

3 


9 

56.57 

6.48 

17 

5 

4 

1 

10 

58.2 

-31.75 

10 

2 

2 


11 

85.08 

-39.58 

10 

2 


2 

12 

43.67 

-85.28 

10 

1 

1 


13 

112.81 

55.61 

8 

2 


1 

14 

138.36 

-10.87 

10 

2 

1 

1 

15 

159.21 

-27.53 

8 

1 

1 


16 

321.97 

-60.77 

16 

3 


3 

17 

316.27 

1.7 

8 

1 



18 

44.13 

-45.61 

8 

1 

1 


19 

60 

-58 

8 

2 

1 

1 

20 

195.95 

28.37 

14 

2 

1 

1 

21 

349.47 

15.55 

16 

4 

3 


22 

197.2 

-18.49 

8 

2 

2 


23 

360 

-31 

8 

3 

2 

1 

24 

17.46 

-45.9 

8 

1 


1 

25 

68.41 

-13.25 

8 

1 

1 


26 

87.51 

-25.09 

8 

1 



Total 




50 

30 

16 


Notes. We group neighboring clusters in non-overlapping analysis regions of 
interest (ROI) that are defined by center location and the associated radius. Note 
that with overlap, we refer to the squared counts map inscribed in the circle 
defined with the coordinates in the table. We list the center of each region along 
with its radius together with its cluster content. Columns from left to right: ROI 
ID, R.A. (J2000), decl. (J2000), radius of ROI, number of total clusters in ROI, 
number of CC-clusters in ROI, and number of NCC-clusters in ROI. 

or annihilation to future studies. For a study of y rays origi- 
nating from star-forming galaxies in galaxy clusters, see, e.g.. 
Storm et al. (2012). Hydrodynamic simulations of galaxy clus- 
ters in a cosmological framework that include CR physics show 
an approximately universal spectral and spatial CR distribu- 
tion within clusters as a result of hierarchical structure growth 
(Pinzke & Pfrommer 2010). The ]/-ray emission induced by 
decaying neutral pions dominates over the inverse Compton 
emission from primary shock-accelerated electrons or secon- 
daries injected in hadronic CR p-p interactions within clusters 
(Miniati 2003; Pfrommer et al. 2008; Pinzke & Pfrommer 2010). 

Using a very simplified analytic model that employs a CR 
power-law energy spectrum dncR/ds oc s~^ with spectral 
index a = 2 (i.e., equal CR energy density per logarithmic 
energy interval), Kushnir & Waxman (2009) claim that the 
IC emission from primary accelerated electrons at accretion 
shocks dominates over pion decay y rays by a factor of 
~150(^e/Cp)(^7'/10keV)“^/^, where and denote the 
fraction of shock-dissipated energy that is deposited in CR 
electrons and protons, respectively. Instead of the centrally 
concentrated pion decay emission, this model would give rise 
to a spatially extended emission reaching out to the accretion 
shocks beyond the virial radius. 

However, there are two simplifying model assumptions that 
conflict with results from numerical cluster simulations and 


observations of non-thermal emission of clusters and super- 
nova remnant shocks, rendering the conclusions about this ap- 
parent dominance of the primary inverse Compton emission 
questionable. 

First, the assumed spectral index is in conflict with numer- 
ical simulations of cosmological structure formation that have 
shown a continuous distribution of Mach numbers with weaker 
flow and merger shocks being more numerous in comparison 
to strong shocks with Mach numbers exceeding M > 6 (Ryu 
et al. 2003; Pfrommer et al. 2006; Skillman et al. 2008; Vazza 
et al. 2009, 2011). This necessarily implies a softer effective 
spectral injection index of primary shock-accelerated particles 
of ofinj —23, which is also consistent with observed spectral 
indices of radio relics. In fact, elongated relics show a mean 
radio spectral index of {a^) = 1.3 (Feretti et al. 2012), which 
translates to a cooling-corrected spectral index of CR electrons 
of Ck: = 2ck:y = 2.6 (which is a lower limit since the uncor- 
rected injection index could be as high as Ckf = 2a^ I = 3.6). 
Hence phenomenologically, the relevant shock strengths for ra- 
dio relic emission are on average characterized by small Mach 
numbers of M < 2 and inconsistent with hard CR indices of 
O' = 2. It can easily be seen that the different power-law indices 
are indeed the reason for the strongly differing conclusions on 
the importance of primary inverse Compton emission. Elec- 
trons with an energy of 500 GeV can Compton upscatter CMB 
photons to ]/-ray energies of around 1 GeV. Adopting the ef- 
fective spectral injection index of primary electrons, ofinj — 2.3, 
and assuming a post-shock temperature at a typical accretion 
shock of kr ~ 0.1 keV, we find a flux ratio of the Kushnir & 
Waxman (2009) model and that by Pinzke & Pfrommer (2010) 
of (500GeV/0.1 keV)°-3 '°-® ~ 800 . . . 6 x 10^. 

Second, ^e/Cp 1 is inconsistent with the observed ra- 
tio for Galactic CRs, which has a differential energy ratio of 
Ke/p — 0.01 at 10 GeV (Schlickeiser 2002) and observations of 
supernova remnants constrain the ratio ^e/?p 0.001 (Edmon 

et al. 2011; Morlino & Caprioli 2012). Assuming universality of 
the shock acceleration process, statistical inferences about the 
CR distributions at the solar circle and non-thermal modeling 
of individual supernova remnants indicate that electron acceler- 
ation efficiencies are very subdominant in comparison to that of 
protons. 

Hence, as a baseline model, we apply the simulation-based 
analytical approach for the CR distribution (Pinzke & Pfrommer 
2010), which only requires the gas density profile inferred 
from X-ray measurements as input. Note that these simulations 
account only for advective CR transport, where the CRs may 
be tied to the cluster plasma via small-scale tangled magnetic 
fields with cored CR profiles as a consequence. Additional 
CR transport such as CR diffusion and streaming have been 
neglected. Furthermore, we neglect the potentially important 
contribution from re- accelerated CRs (e.g., Brunetti & Lazarian 
2007, 2011). 

The pion decay y-ray flux above energy £" as a result of 
hadronic CR interactions, F^g^p(>F^), can be parameterized as 

= A, (>E) J^dVK„.^y(R), (3) 

where the integral extends over the cluster volume V, 
kj^o^y (> E) is the spectral y-ray distribution, and Kj^o^y(R) 
is the spatial distribution, both given in Pinzke & 
Pfrommer (2010). The parameter, , is a dimensionless univer- 

sal scale factor, common to all the clusters in the (sub)sample. 
The predicted y-ray flux above the minimum energy threshold 
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Figure 3. Expected surface brightness profiles for the three spatial models 
considered in this analysis. We show the prohles for two clusters, Coma (black 
lines) and A400 (red lines), which have comparable distances (z = 0.02) such 
that the flux difference corresponds to the difference in mass. We note that 
the ICM model (dashed line) only shows small differences with respect to our 
simulation-based baseline model (solid line). In contrast, the flat CR pressure 
model (dashed-dotted line) implies a flattening toward the outskirts of the cluster. 
(A color version of this figure is available in the online journal.) 


500 MeV, > 500 MeV) is calculated using the formal- 

ism in Pinzke et al. (2011). We tabulate these values in Table 1. 
The quoted values of > 500 MeV) correspond to a 

maximal efficiency ^p,max = 0.5 for diffusive shock accelera- 
tion of CR ions at structure formation shocks which translate 
into Ay = I with correspondingly smaller values of for 
smaller efficiencies (obeying however a nonlinear relation). For 
completeness, we show the CR formalism in Appendix B. 

The cluster brightness profile is used to fit the emission from 
each cluster and is derived from the line-of-sight (l.o.s.) integral 
of the ]/-ray emissivity 

> E) = Ay XjrO^y (>E) 



where I is the l.o.s. distance in the direction \l/ that the 
detector is pointing and R(l) = — 2DJ cos 'F is 

the cluster radius. Here Da is the angular diameter distance 
from the Earth to the center of the cluster halo and cos 'F = 
cos 0 cos xj/ — cos (p sin 0 sin i/^, with 0 being the azimuthal angle 
and 0 the polar angle. The angular integration dQ. = sin OdO dcp 
is performed over a cone centered around 0 . The spatial features 
of our model are described in more detail in Appendix B. 

Outside the very center (r > 0.03^200)5 this model predicts a 
rather fiat CR-to-thermal pressure profile, i.e., (Xcr) const. 
Most of the emission is contributed from the region around 
the core radius, which is well outside the central parts. In 
order to compare the chosen spatial CR profile, we contrast 
our analysis of the simulation-based model with two additional 
configurations in which the CR profile is derived from a constant 
XcR profile (ICM model) and a constant Pcr profile (flat model). 
The normalizations of these CR profiles are fixed by assuming 
that the total CR number within R200 in our baseline model 



Mi5 = M2oo/10*'^M® 


^ 0.025 r 

S 0.020 h 

d 0.015 h 
^ 0.010 r 



R! R200 


Figure 4. Relative cosmic-ray pressure Xcr within radius R. We show the 
ratio between cosmic-ray pressure and thermal pressure for 14 simulated galaxy 
clusters with different mass. In the upper panel, the color scheme shows the 
relative pressure within different radii, in descending order from 1.0 x R 200 
to 0.2 X R 200 (equally spaced in log R). Each simulated cluster is denoted by 
a X and the mass dependence of Xqr as a function of radius are denoted by the 
solid lines. The middle panel shows the radial dependence of the normalization 
of Xqr. The lower panel shows the radial dependence of the slope of the mass 
dependence of Xcr. We And that the relative pressure increases as a function of 
the radius, but it decreases with increasing cluster mass. 

(A color version of this flgure is available in the online journal.) 

is conserved.^^ For illustration purposes, we show the surface 
brightness profiles using our three spatial emission models in 
Figure 3 for the case of the (massive) Coma cluster and the much 
less massive cluster A400. 

In our framework, the derivation of A^ also allows us to con- 
strain the CR-to-thermal pressure ratio, Xqr, in the ICM using 
the virial mass and virial radius of each cluster in our sample, 
since Ay a (Xcr), where (Xcr) = (Pcr) v / (Pth) v and the 
brackets indicate volume averages. To this end, we make use 
of the set of 14 galaxy cluster simulations presented in Pinzke 
& Pfrommer (2010), which span almost two orders of mag- 
nitude in mass and include different dynamical states ranging 
from relaxed to merging clusters. We show the CR-to-thermal 
pressure ratio Xqr as a function of radius and cluster mass in 
Figure 4. Xqr decreases for smaller radius approximately in- 
versely with gas temperature since a composite of CRs and 
thermal gas favors the gas component over CRs upon adia- 
batic compression (e.g., Pfrommer & EnBlin 2004a). At a fixed 


Note that CR streaming and diffusion — as spatial transport 
processes — conserve the total number of CRs. However, outward streaming 
changes their number density as a function of radius and transforms a peaked 
CR profile into an asymptotically flat one. Hence, to map an 
advection-dominated profile (i.e., our simulation-based baseline model) to the 
corresponding asymptotically flat profile that results from streaming, we 
compute the volume integral of the number density before and after CR 
streaming and normalize the latter such that the total number of CRs is 
conserved. 
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radius, Xqr has a negative trend with mass, that is mainly driven 
by the virial temperature scaling of the thermal pressure dis- 
tribution. We have (Xcr) oc (C)/(Pth) C/kT2oo ~ , 

where C is the normalization of the CR distribution function and 
C = Cnipl p denotes the dimensionless normalization, which 
scales with cluster mass as (C)y a (Pinzke & Pfrommer 
2010) and partially offsets the virial mass scaling of the cluster 

2/3 

temperature k72oo oc 

To formalize these considerations, we fit an empirical relation 
of Xcr{R/ R im, ^ 20 o) to the simulated data points and obtain 


Xcr = 0.023 Ay 


R M200 

R200) Vio'^Mq 


-o-239(45)- 


(5) 


where Ay refers to our universal scale factor derived in the joint 
likelihood analysis. We tabulate the values for (Xcr) for two 
different integration radii, R 200 and the half-light radius Rhl 
assuming = 1 in Table 1.^^ 


4.2. Background Model 

The background model for each ROI in this analysis includes 
templates for the diffuse Galactic and isotropic emission com- 
ponents as well as individual y-ray sources reported in the 2nd 
Fermi-LKV catalog (Nolan et al. 2012; 2FGL). We model ex- 
tended 2FGL sources according to the spatial templates provided 
by the Fermi Science Support Center. Unless stated otherwise, 
we use the standard diffuse and extragalactic y-ray background 
templates that are recommended for performing data analysis of 
reprocessed LAT data.^^ We note that the Galactic diffuse emis- 
sion model we use includes a residual component of diffuse 
]/-ray emission that is not modeled by any template. This com- 
ponent is smoothly varying and does not contribute importantly 
to the intensity >500MeV in the regions considered here.^^ 

In the background model for each ROI, we include the union 
of 2FGL sources within the ROI radius enlarged by 5° and 2FGL 
sources located within 10° of any galaxy cluster in the ROI. 



ROI 


Figure 5. Fit results for the normalization for the Galactic diffuse emission 
(top) and the isotropic diffuse emission (bottom) in each of the 26 analyzed 
ROIs. The dot-dashed lines indicate a nominal value of one associated to diffuse 
templates exactly modeling the emission in the regions. Our regions show a 
narrow scatter for the Galactic diffuse emission and a slightly larger scatter for 
the isotropic diffuse component. The latter is also associated with a minor bias 
toward normalization values > 1 . 

likelihood function introduced in Section 3.1 in order to avoid 
an overflow of free parameters and to ensure convergence of the 
maximization. The normalizations of the Galactic and isotropic 
diffuse components are left free in each ROI for the joint 
likelihood fit. 


5. RESULTS AND DISCUSSION 

We perform our analysis first by treating all 50 clusters in one 
common set. We call this the combined sample. Recalling the 
discussion in Section 2, we then separately investigate CC and 
NCC clusters. 


4.3. Free Parameters of the Background Model 

Since the average virial radius of clusters in our sample is 
less than 2°, we leave the normalizations of all sources free 
that are contained within a 4° radius around each cluster. In 
addition, we allow the normalization of the templates used to 
model the Galactic foreground and isotropic diffuse emission to 
vary freely. 

One shortcoming of using the 2FGL catalog (based on 2 yr 
of LAT observations) to search within a data set covering 4 yr is 
that spectral parameters (in particular for variable sources) may 
have substantially changed. To account for this variability, we 
free the normalization of sources that coincide with bright spatial 
residuals, and we determine their values through performing a 
fit using a background-only model to obtain the best-fit for the 
null-hypothesis. 

This procedure produces a large number of free parameters 
which are then fixed to their maximum likelihood values from 
the background-only model fit when maximizing the joint 


Here we define the half-light radius as the radius within which half of the 
emission originates. Note that this radius is usually smaller than R 200 / 2 . 

For the Galactic diffuse emission, we use the template gll_iem_v05 and 
for the isotropic y-ray background iso_clean_v05. 

Further details on this new model can be found on the FSSC Web site. 


5.1. Background-only Fit 

Our ROIs are well described individually by the null hypoth- 
esis, i.e., despite the increase in data volume, our results are 
consistent with emission from only previously detected indi- 
vidual Fermi-hAI sources and diffuse emission provided by 
the Galactic and extragalactic diffuse models, respectively, with 
the exception of three ROIs: the ROIs containing A400 and, 
less prominently A1367 and A3 112 exhibit residual emission 
located within the virial radius of each respective cluster. We 
leave these excesses unmodeled in our baseline analysis and ad- 
dress the interpretation of these residuals in Section 5.3. It is also 
reassuring that the fitted normalizations of the two global diffuse 
backgrounds are narrowly scattered around the nominal value of 
1 for both the extragalactic (isotropic) and the Galactic diffuse 
component (refer to Figure 5) across our entire sample. The 
isotropic component shows a slight bias toward normalizations 
> 1 which however has negligible effect on our results. Figure 6 
shows a stacked residual significance map for the full sample and 
CC and NCC sub-samples. These residual significance maps are 
created from theoretical model maps of predicted counts from 
the best-fit null hypothesis model summed over each cluster lo- 
cation. The combined model maps M are then subtracted from 
the stacked counts maps C from the same region and the resid- 
ual significance R is computed sls R = (C — M)/\fM. For the 
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Figure 6. Spatial residual significance maps for the combined sample (top) and 
the two sub-samples (bottom). The dashed circle corresponds to a radius of 1?5 
which covers the majority of the emission region, assuming that the emission is 
contained within the cluster’s virial radius. The stacked maps are created using 
the galaxy cluster centers listed in Table 1 . We do not observe any significant 
excess in the combined sample. If observed, any excess in these maps would 
be extended on at least the scale of the effective point-spread function (PSF) of 
LAT. Each pixel has a size of 0?25. 


spectral residuals, we refer the reader to Appendix E. No obvi- 
ous excess is visible in either the spectral or spatial residuals. 
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Figure 7. Distribution of ROI TS values for full sample (blue solid line), cool 
core (red dashed line), and non-cool core (green dash-dotted line). We discuss 
the notable exceptions with TS ^ 9 in Section 5.3. 

(A color version of this hgure is available in the online journal.) 


Table 3 

Joint Scale Factor Limits 


Model 

Combined 

CC 

NCC 

Pinzke & Pfrommer (2010) 

0.40 (7.4) 

0.49 (4.7) 

0.44 (2.4) 

Constant Xcr (ICM model) 

0.48 (7.2) 

0.64 (4.9) 

0.49 (2.4) 

Constant PcR (flat model) 

1.78 (9.7) 

3.02 (5.2) 

1.71 (5.0) 


Notes. Summary of joint scale factors. Columns 2-4 indicate the 95% upper 
limit on the joint scale factor Ay for the respective samples. The global TS 
values of each setup and sample are given in brackets. 


5.2. Global Significance and Constraints 
on Common Scale Factor Ay 

We then repeat the fitting procedure including a model of 
our galaxy clusters with the predicted y-ray flux > 

500 MeV) and using the spatial template and spectral form 
proposed by Pinzke et al. (2011), leaving only to vary 
freely. We show the distribution of associated TS values for 
the respective samples in Figure 7. 

Assuming that backgrounds are properly modeled, and that 
CR physics governing the ]/-ray emission of clusters is indeed 
universal, we calculate the best-fit value of the combined scale 
factor for the full sample along with the two morphological 
sub- samples. The global TS value of the scale factor for the 
full sample of 50 clusters is 7.3, corresponding to a 2.1a 
evidence. We note that the largest contributors to this signal 
are however the aforementioned excesses which are spatially 
coincident with A400 but also less prominently from excesses 
toward A1367 and A31 12. We discuss these special cases in the 
next section. Removing these clusters from the sample results 
in a drop of the significance below 2a, yielding an upper limit 
to the common scale factor = 0.29 at 95% C.L. for 
the whole cluster sample containing the remaining 47 galaxy 
clusters. While individually the excess toward A400 yields a 
higher significance than the excesses toward either A3 112 or 
A1367, in the combined limit, the contribution of this excess 
becomes negligible due to the lower flux prediction (which 
mainly determines the weight assigned to each cluster in the 


joint analysis), as compared to, e.g., the contribution from Coma 
that dominates the upper limit. 

Since at present we can neither claim nor refute the origin of 
the observed excesses being due to y rays from the ICM, we 
calculate upper limits on the universal scaling factor, leaving 
those respective excesses unmodeled. For our whole sample, we 
find a combined limit of A^^ = 0.41 at 95% C.F. Considering 
the NCC/CC subsamples, we find weaker limits A^^ = 0.47 on 
the scale factor for NCC systems as compared to Aj/^ = 0.49 
for the CC subsample. 

We also calculated the limit on the combined scale factor for 
our two alternative spatial CR profiles. For the ICM model, we 
obtain Aj/^ = 0.48, and for the fiat model, the combined limit 
is roughly a factor of four larger with respect to the results from 
the baseline model, yielding A^^ = 1.78, at 95% C.L. The 
associated global TS values are 7.2 and 9.7, respectively.^^ In 
addition, a fiat CR profile is preferred in the case of the Coma 
cluster which is also the cluster that contributes most to the 
constraints in the NCC cluster subsample. We provide the final 
values for our three setups in Table 3. 

For the combined sample, we exclude Ay = 1 at more than a 
5 a confidence level, while for CC/NCC, it is excluded at more 
than a 4 a confidence level. 


Removing A400, A1367, and A31 12 from the cluster sample yields smaller 
global TS values of 2.8 for the ICM model and 4.7 for the flat model, 
respectively. 
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Figure 8. TS maps from an unbinned search in a 5° x 5° region centered on each of our notable clusters, A1367, A3 112, and A400. All excesses are found within 
the assumed cluster virial radius (dashed white circle), albeit marginally offset from the respective cluster centers (0?3 for both A400 and A1367 and 0? 1 for A31 12). 
Each pixel has a width of 0? 1 . The white x indicates the best-fit position of a previously detected 2FGL point source. 

(A color version of this hgure is available in the online journal.) 


Table 4 

Best-fit Positions of Excess Emission 


Host Cluster 

R.A. 

(°) 

Deck 

(°) 

^68 

(°) 

A0400 

44.68 

5.86 

0.03 

A1367 

176.25 

19.54 

0.04 

A3112 

49.56 

-44.22 

0.02 


Notes. We report the best-fit positions from our refined search 
using the gttsmap tool that employs a maximum likelihood 
analysis in order to localize a new point source. The columns 
from left to right: name of the host cluster, R.A. (J2000), 
deck (J2000), and uncertainty All values are given in 
degrees. 


5.3. The Case ofA1367, A3 112, andA400 

We obtained a TS value of 13.8, corresponding to a pre- 
trial factor significance of 3.1a individually for each of the 
candidate y-ray sources at the locations of A1367 (Leo cluster) 
and A3 112. Conservatively, assuming a binomial distribution 
for the trial probability, we find that these excesses correspond 
to a post-trial significance of 2.6a. For each of these regions, we 
calculated a TS map using an unbinned likelihood method in a 
5° X 5° region centered around each cluster with a grid spacing 
of 0? 1 between test positions. The TS value of a putative point- 
like source with a hadronic CR spectrum is evaluated at each 
test position on the grid to create a spatial map of the excess 
emission. We show these TS maps in Figure 8. 

In both cases, we find that the excess emission, albeit 
marginally offset from the center of the cluster (0?3 for A 1367 
and 0?1 for A3 11 2) is still contained within the virial radius of 
the respective cluster. For A 1367 the difference in TS obtained 
at the center of the cluster and the peak TS position as shown 
in Figure 8 is 15, while for A31 12 the offset is smaller than the 


resolution used to create the TS maps. In order to test whether the 
emission is more appropriately modeled assuming an extended 
emission profile over a new point source, we follow the analysis 
presented in Lande et al. (2012), which gauges the spatial extent 
of the source by taking the difference, TSgxt between the TS 
for a point source signal hypothesis and the extended source 
signal hypothesis. Lande et al. (2012) found that sources with 
TSext > 16 are confidently ascertained to be spatially extended 
beyond the LAT point-spread function (PSF).^^ For A1367 and 
A3 112, we find TSext = 2.6 and TSext = 0-9, respectively. 
Assuming that the excesses toward all three clusters constitutes 
point-like emission, we first obtained a better estimate for the 
location of the excess using the gtfindsrc tool and then 
repeated the calculation of TS maps using a finer binning of 
0?02. We report the best-fit values of these excesses along with 
their uncertainties in Table 4. 

We also investigate the spectral behavior of the excesses 
toward A 1367 and A3 112, by replacing the hadronic CR 
spectrum with a featureless power-law, such that the flux 
becomes: 


fzPL 


(£) = A, X F“p(>500MeV) 


(1 - n X E~^ 


( 6 ) 


where F is the spectral index and L^min = 500 MeV and 
= 200 GeV correspond to the energy range of the 
analysis. > 500 MeV) denotes the expected integral 

flux above ^min- correspouds to the scale factor introduced 
in Section 5.2. For this test, we leave both Ay and F free to vary. 
We report the best-fit values of these parameters in Table 5. 

For both A1367 and A3 112, we find that harder spectral 
indices are preferred over softer, with a best-fit value for 

In Lande et al. (2012), the authors used a disk to perform their studies. 


Table 5 

Spectral Model Comparison 


Cluster 

Ay 

(Hadronic Model) 

TS 

(Hadronic Model) 

Ay 

(Power Law) 

r 

(Power Law) 

TSpL 

(Power Law) 

A1367 

3-1 

13.8 

1.7 ±0.9 

1.7 ±0.3 

17.3 

A3112 

3±2 

13.8 

2.0 ± 1.5 

1.7 ±0.4 

16.1 

A400 

OQ+ll 

52.7 

43 ±8 

2.3 ±0.2 

52.8 


Notes. Spectral model comparison of clusters that exhibit excess emission. Shown are the best-fit values for Ay along with their 
associated TS values for the hadronic model and the corresponding values for Ay when replacing the spectrum by a featureless 
power-law of index E, given by Equation (6). The last column indicates the obtained TS value with the power-law fit. 
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r = 1.7 zb 0.3 and T = 1.7 zb 0.4, for A1367 and A3112, 
respectively. 

While none of the aforementioned tests decisively exclude 
the attribution of the observed excesses to CR-induced y-ray 
emission from the ICM, we note that both A3 112 and A1367 
are hosts to head-tail radio sources which may be the source of 
the observed ]/-ray emission (see, e.g., Gavazzi & Jaffe 1982, 
1987 for A1367 and Costa & Loyola (1998) for A3112). A 
discussion of supporting multifrequency arguments is given in 
Appendix D. 

The excess found in A400 yields a TS value of 52.7 which 
nominally corresponds to a significance of 7.3a pre-trial (6.7a 
post-trial). We performed the same tests as in the case of A1367 
and A3 1 12 and find that the excess is contained within the cluster 
virial radius, although the TS map indicates that the emission is 
about 0?3 offset from the cluster center. The difference between 
the TS evaluated at the cluster center and at the position of 
the excess is 38. The test for extendedness (centered at the 
cluster center) yields TSext = 15.0, indicating a preference for 
an extended source, which is likely explained by the offset of the 
excess. This casts further doubt on the explanation of the excess 
being ICM emission that moreover would be concentrated 
toward the cluster center. Fitting the excess with a power-law 
spectral model as in the case of A3 112 and A1367 yields a 
best-fit value for F = 2.3 zb 0.2 and Ay = 43 zb 8 with an 
associated TS value of 52.8, which is similar to that obtained 
for the hadronic model. 

In addition to these tests, we searched for source variability 
using aperture photometry and found no indications of variabil- 
ity on time scales of one month. We note, however, that the 
obtained scale factor for A400, Ay = 39t.\o’ strong tension 
with baseline model expectations and the scale factor constraints 
derived from other clusters in our sample. If the excess toward 
A400 constitutes a signal, the calculated upper limit from A400 
corresponds to the usual statement that scale factors larger than 
the upper limit are inconsistent with the data on the stated con- 
fidence level. If the excess instead stems from the unmodeled 
background, then the upper limit means that scale factors larger 
than the upper limit would make the model more inconsistent 
with the data than the background-only hypothesis allows at the 
stated confidence level. In both cases, the result is a conservative 
upper limit and in addition (as mentioned in Section 5.2) A400 
does not affect the combined upper limit sizeably. Removing 
A400 from the sample yields a marginally stronger upper limit 
on the common scale factor for the combined sample of the 
remaining 49 clusters of A^^ = 0.40. 

In the conservative CR model of Pinzke & Pfrommer (2010), 
Ay = 1 corresponds to a pion decay y-ray flux owing to 
CR protons accelerated at structure formation shocks with 
a maximal acceleration efficiency and neglecting active CR 
transport.^ ^ Allowing for CR streaming transport would cause a 
net CR flux to the dilute outer cluster regions and would reduce 
the }/-ray yield for that acceleration efficiency. Additionally, 
Pinzke & Pfrommer (2010) presented an optimistic CR model 
including effects which enhance the predicted ]/-ray yield by a 
factor of 2-3 depending on the cluster mass.^^ However, that 
model does not account for CRs that are injected by AGNs 

The acceleration efficiency that was used in the simulations was based on 
observations of supernova remnants and theoretical calculations of diffusive 
shock acceleration. It is unlikely that there are more efficient mechanisms at 
work. 

A part of the enhancement was due to the numerical limitation of smoothed 
particle hydrodynamics to excite Kelvin-Helmholtz instabilities and fully mix 
a ram-pressure stripped interstellar medium into the ICM. 


over the cluster lifetime, which could also produce pions in 
inelastic collision with the ICM. The total energy dissipation 
by gravitational shocks exceeds that of AGNs for large clusters, 
making it unlikely that AGN-injected CRs dominate the diffuse 
}/-ray signal. However, this argument does not apply for smaller 
CCs (in particular in their central cluster regions) where the 
AGN appears to dominate the energy budget and could possibly 
also give rise to observable y-ray emission (see, e.g., Pfrommer 
2013, for the interpretation of the low-state of M87 in terms of 
diffuse pion decay emission while the high state is attributed 
to jet-induced emission). Hence, the most likely explanation for 
the possible signal with Ay ~ 39 in A400 is jet-related emission 
from point sources projected onto or within the cluster. A 
potential source for the emission in A400 may be the quadruple 
head-tail system 3C 75 which also shows X-ray core emission 
in the galaxies (Owen et al. 1985; Hudson et al. 2006). However, 
given that 3C 75 is located toward the center of A400 and the 
excess is 0?3 offset, which is about eight times larger than the 
68% error radius for a pointsource, this possibility is unlikely. 

5.4. Individual Upper Limits on the y-ray Flux 

Assuming that each cluster in our sample can be modeled 
according to our description in Section 4.1, one can also derive 
individual limits on / where i refers to the individual cluster. 
These individual limits can be propagated into flux limits. 

In addition, in order to assess the impact of modeling the 
clusters in our sample as extended sources, we have derived 
individual flux upper limits modeling the clusters as point 
sources. 

In Figure 9, we show these two cases, contrasting the indi- 
vidually derived y-ray flux limits for extended emission from 
those derived when assuming the cluster emission to be point- 
like. We tabulate the former for various energy bands along with 
their associated (pre-trial) TS values in Table 6. We note that the 
limit on A^^ , derived from the Coma cluster alone is comparable 
to the jointly derived limit from the full sample of all 50 clusters 
in this study, emphasizing its weight in the joint analysis. We 
investigated the potential dependence of the upper limits on the 
Galactic diffuse emission model. In the great majority of cases 
the limits vary by less than 30% for the range of models that 
we considered. Those clusters for which the dependence on the 
model was more sensitive are marked in Table 6. However, this 
sensitivity does not affect the combined limit. 

In Figure 10, we show individual y-ray flux upper limits 
that are derived for the CR profile following: (1) a constant 
XcR profile (ICM model) and (2) a constant Pcr profile (fiat 
model). We tabulate these values along with their (pre-trial) 
TS values in Table 7. We find that the limits derived from 
the first model are very similar to those with our simulation- 
based model. On the other hand, the fiat model (i.e., constant 
CR pressure profile) yields less stringent upper limits. We also 
note that the associated TS values for the fiat CR profile are 
generally marginally higher than either those derived from the 
thermal ICM or our simulation-based approach. This is expected 
since the choice of a fiat CR profile yields a substantially flatter 
y-ray surface brightness profile, which in turn provides a better 
fit to the data compared to a cored profile, in particular, when 
considering that the excess emission we report in the previous 
section is offset from the cluster center. We also note that from 
the three excesses found, only those toward A3 112 and A400 
remain with TS > 9 when using the fiat CR profile instead. 

The flux limits we derive substantially improve over previous 
limits (Ackermann et al. 2010b) due to the increase in data 


12 


The Astrophysical Journal, 787:18 (26pp), 2014 May 20 


Ackermann et al. 




Figure 9. Shown are the 95% upper limits on hadronic CR-induced y-ray flux for each of our 50 galaxy clusters in this analysis. We show the individually derived 
upper limits for both the extended emission (red downward triangle) and assuming the cluster emission to be point-like (blue circle). 

(A color version of this flgure is available in the online journal.) 


volume and improved modeling of the y-ray sky as well as 
improved instrument understanding reflected by the use of 
reprocessed LAT data with on-orbit calibrations. Finally, we 
note that in particular for spatially extended clusters such as 
Coma or Fornax, the limits are substantially weakened with 
respect to when modeling them as point-like objects. 

Motivated by the Fermi-LP^ detection of few bright clus- 
ter galaxies (e.g., 2FGL J0627. 1—3528 in A3392 and 2FGL 
J1958.4-3012 in RXC J1958— 3011), a recent stacking study 
by Dutson et al. (2013) investigated a large sample of galaxy 
clusters that was selected according to the radio flux of bright 
cluster galaxies. Although based on a different scientific prior 
and methodology than our cluster analysis, the determined flux 
limits can be compared to the point- source upper limits re- 
ported here. About a dozen clusters are in common between both 
studies. However, Dutson et al. (2013) used the respective co- 
ordinates of the bright cluster galaxy, which are not necessar- 
ily consistent with the cluster center coordinates considered 
in our studies. Given the LAT PSF (Bregeon et al. 2013) and 
the considered ROIs this does not constitute a severe handicap 
for comparison. The use of different exposures (45 months in 
Dutson et al. 2013 and 48 months in this work), the use of differ- 
ent source models, and, perhaps most particularly noteworthy, 
the use of reprocessed LAT data with associated different Galac- 
tic and isotropic diffuse models (gal_2yearp7v6_v0 versus 
gll_iem_v05 and iso_p7v6source versus iso_clean_v05 


respectively), as well as different analysis energy thresholds, 
render a strict comparison more problematic. For the major- 
ity of common clusters, the limits in Dutson et al. (2013) are 
marginally less sensitive, as expected regarding the slightly less 
exposure and the rather moderate changes in the diffuse back- 
ground models. 

However, there are two noticeable exceptions: A85 and 
A2634 appear to have more constraining upper limits in Dutson 
et al. (2013), besides less exposure and a lower analysis 
threshold. The discrepancies could be explained by differences 
in the construction of the ROI (treatment of variable sources, 
sources to be too faint to be in the 2FGL catalog, size of ROI). 
Taking the respective flux limits at face value, the differences do 
not amount to more than 35% between both studies. Our limits 
on extended cluster emission cannot be meaningfully compared 
to the point-source upper limits in Dutson et al. (2013) as they 
constitute alternative scientific priors for a different scientific 
problem. However, our individual limits on the ]/-ray flux, while 
being specifically derived within the framework of the universal 
CR model by Pinzke & Pfrommer (2010), can, in principle, be 
used to constrain other classes of models. 

5.5. CR-to-Thermal Pressure Ratio (Xcr) 

We show the resulting upper limits on the CR-to-thermal pres- 
sure ratio, {Xcr) in Figure 11. These numbers were obtained 
by scaling the {Xcr) values in Table 1 with the limit on 
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Cluster 

aUL 

UUL 

^y,500MeV 
(x 10-10) 

UUL 

^7,1 GeV 

(xlO-11) 

uUL 

^y,10GeV 
(x 10-12) 

T8 

aUL 

uUL 

^y,500MeV 
(x 10-10) 

UUL 
^7,1 GeV 

(xlO-^1) 

UUL 

^7,10GeV 

(xlO-12) 

TS 


Extended 

Extended 

Extended 

Extended 

Extended 

Point-like 

Point-like 

Point-like 

Point-like 

Point-like 

2A0335 

0.52 

1.5 

6.7 

3.6 

0.0 

0.45 

1.3 

5.9 

3.1 

0.0 

A0085 

1.35 

3.9 

18.0 

9.5 

1.0 

1.26 

3.7 

16.9 

9.0 

1.3 

A0119 

3.54 

5.0 

23.0 

12.2 

2.7 

2.65 

3.8 

17.2 

9.1 

0.8 

A0133 

2.32 

1.7 

7.6 

4.0 

0.0 

2.23 

1.6 

7.3 

3.9 

0.0 

A0262 

1.67 

2.0 

9.3 

4.9 

0.0 

1.07 

1.3 

5.9 

3.1 

0.0 

A0400 

50.11 

22.2 

101.7 

53.6 

52.7 

41.20 

18.3 

83.6 

44.1 

37.7 

A0478'' 

1.14 

2.8 

12.7 

6.7 

0.0 

0.88 

2.1 

9.8 

5.2 

0.0 

A0496^ 

1.94 

5.5 

25.2 

13.3 

1.9 

1.46 

4.1 

18.9 

10.0 

0.5 

A0548e 

10.78 

2.7 

12.3 

6.5 

0.1 

9.02 

2.2 

10.3 

5.4 

0.0 

A0576 

4.39 

2.7 

12.6 

6.6 

0.2 

3.80 

2.4 

10.9 

5.7 

0.3 

A0754 

2.86 

4.8 

22.2 

11.7 

1.1 

2.49 

4.2 

19.3 

10.2 

0.7 

A1060 

1.89 

4.3 

19.7 

10.4 

0.4 

1.28 

2.9 

13.3 

7.0 

0.1 

A1367 

4.08 

7.8 

35.9 

19.0 

13.8 

3.31 

6.4 

29.1 

15.4 

11.2 

A1644 

2.69 

3.5 

16.0 

8.5 

0.2 

2.41 

3.1 

14.3 

7.6 

0.4 

A1795^ 

0.42 

1.3 

5.8 

3.1 

0.0 

0.42 

1.3 

5.7 

3.0 

0.0 

A2065 

5.00 

4.6 

21.1 

11.2 

2.3 

4.69 

4.3 

19.8 

10.5 

2.2 

A2142 

0.47 

1.6 

7.4 

3.9 

0.0 

0.46 

1.6 

7.3 

3.9 

0.0 

A2199 

1.45 

4.3 

19.8 

10.5 

1.9 

1.32 

4.0 

18.1 

9.6 

1.7 

A2244 

1.74 

1.3 

6.0 

3.2 

0.0 

1.63 

1.2 

5.6 

3.0 

0.0 

A2255 

5.21 

4.4 

20.2 

10.6 

5.2 

4.99 

4.2 

19.3 

10.2 

6.3 

A2256 

0.50 

1.3 

5.8 

3.1 

0.0 

0.42 

1.1 

4.9 

2.6 

0.0 

A2589^ 

3.96 

2.7 

12.2 

6.5 

0.0 

3.36 

2.3 

10.4 

5.5 

0.0 

A2597 

1.27 

1.0 

4.4 

2.3 

0.0 

1.06 

0.8 

3.7 

2.0 

0.0 

A2634 

5.95 

3.7 

17.0 

8.9 

0.3 

4.16 

2.6 

11.9 

6.2 

0.0 

A265T 

2.34 

1.9 

8.9 

4.7 

0.0 

1.75 

1.4 

6.6 

3.5 

0.0 

A2734 

2.71 

1.2 

5.5 

2.9 

0.0 

2.58 

1.1 

5.2 

2.7 

0.0 

A2877 

1.87 

0.9 

4.3 

2.3 

0.0 

1.55 

0.8 

3.6 

1.9 

0.0 

A3112 

5.37 

6.0 

27.3 

14.4 

13.8 

5.06 

5.6 

25.7 

13.6 

12.8 

A3158 

1.26 

1.5 

6.9 

3.7 

0.0 

1.17 

1.4 

6.4 

3.4 

0.0 

A3266 

1.24 

3.9 

17.7 

9.4 

2.2 

1.15 

3.6 

16.4 

8.7 

3.6 

A3376 

6.20 

4.3 

19.5 

10.3 

0.0 

3.59 

2.5 

11.3 

6.0 

0.0 

A3822 

4.77 

2.1 

9.8 

5.2 

0.0 

4.22 

1.9 

8.7 

4.6 

0.0 

A3827 

0.66 

0.6 

2.9 

1.5 

0.0 

0.61 

0.6 

2.7 

1.4 

0.0 

A392U 

4.63 

2.2 

9.9 

5.2 

0.1 

4.11 

1.9 

00 

bo 

4.6 

0.0 

A4038 

1.84 

2.5 

11.3 

6.0 

0.1 

1.48 

2.0 

9.1 

4.8 

0.0 

A4059 

2.04 

2.0 

9.1 

4.8 

0.0 

1.86 

1.8 

8.2 

4.4 

0.0 

COMA^ 

0.35 

4.0 

18.5 

9.8 

0.7 

0.22 

2.5 

11.3 

6.0 

0.1 

EXO0422 

0.70 

0.5 

2.3 

1.2 

0.0 

0.65 

0.5 

2.2 

1.1 

0.0 

FORNAX 

3.73 

3.1 

14.3 

7.6 

0.0 

1.29 

1.1 

4.9 

2.6 

0.0 

HCG94 

6.34 

2.4 

11.0 

5.8 

0.0 

5.17 

2.0 

9.0 

4.7 

0.0 

HYDRA-A 

2.57 

4.3 

19.6 

10.4 

3.1 

2.44 

4.1 

18.6 

9.8 

3.2 

IIIZw54 

10.97 

5.0 

22.7 

12.0 

0.5 

9.74 

4.4 

20.2 

10.6 

0.3 

IIZwl08 

3.31 

1.3 

6.1 

3.2 

0.0 

3.25 

1.3 

6.0 

3.2 

0.0 

NGC 1550 

3.37 

2.3 

10.3 

5.5 

0.0 

2.70 

1.8 

8.3 

4.4 

0.0 

NGC 5044 

4.96 

3.0 

13.8 

7.3 

0.0 

4.29 

2.6 

11.9 

6.3 

0.0 

RXJ2344 

5.07 

2.9 

13.3 

7.0 

1.3 

5.10 

2.9 

13.4 

7.1 

2.3 

8405^ 

3.11 

1.2 

5.7 

3.0 

0.0 

2.36 

0.9 

4.3 

2.3 

0.0 

8540 

11.47 

4.0 

18.5 

9.7 

1.2 

10.35 

3.6 

16.7 

00 

bo 

0.8 

UGC03957 

8.00 

4.0 

18.2 

9.6 

0.4 

7.31 

3.6 

16.6 

00 

bo 

0.3 

ZwC11742 

1.94 

2.3 

10.4 

5.5 

0.0 

1.77 

2.1 

9.5 

5.0 

0.0 


Notes. The columns contain from left to right the 95% upper limit on Ayj for each cluster, the derived flux upper limit above (500 MeV, 1 GeV, and 10 GeV), the 
associated TS value as well as the same quantities assuming the clusters to be modeled by a point source at the cluster position as given in Table 1 . Fluxes are given in 
photon s“^ cm“^. 

^ The upper limits on Ayj derived using our alternative diffuse models (Section 6.3) for these clusters varied by more than 30% relative to those obtained using the 
standard diffuse emission model. 


from Section 5.2. This procedure assumes universality of the 
CR distribution as suggested by hydrodynamical cosmological 
simulations of clusters (Pinzke & Pfrommer 2010) and implic- 
itly asserts that the active CR transport does not appreciably 
modify the spatial distribution of the CRs. This is justified since 
the impact of CR streaming on the CR distribution of a cosmo- 


logical cluster is not clear to date. Depending on the microscopic 
plasma physics that sets the CR streaming speed (i.e., compet- 
ing damping mechanisms of the CR Alfven waves) and the 
macroscopic distribution of cluster magnetic fields (Pfrommer 
& Dursi 2010; Ruszkowski et al. 2011, for evidence of radial 
bias of the magnetic geometry), CR streaming could either be a 
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Figure 10. Same as Figure 9 but for spatial CR profiles following a constant Acr profile (ICM model, red squares) and a constant PcR profile (fiat model, green stars). 
To allow for an easier comparison, we show the limits from the baseline analysis (Figure 9) in horizontal gray lines. 

(A color version of this figure is available in the online journal.) 


perturbation to the peaked advection-dominated CR distribution 
or cause a substantial flattening by a substantial net outward CR 
flux (EnBlin et al. 2011; Wiener et al. 2013). 

The median upper limit on the CR pressure ratio is (Xcr) < 
0.006 for the combined sample and the NCC subsample within 
Rhl- Those constraints are relaxed to (Xcr) < 0.012 (combined 
sample) and (Xcr) < 0.013 (NCC) within R 200 - Fof CC 
clusters, this limit is less stringent, yielding {Xqr) < 0.008 
and (Xcr) < 0.014 within Rul and R 200 , respectively. 

These limits are more constraining than the previous limits 
on the CR pressure that were obtained through flux upper lim- 
its on individual objects; in particular, they are constraining for 
individual limits on clusters using the initial 18 months of LAT 
data (Ackermann et al. 2010b); also, the limits improve those 
constraints that use four years of Fermi data on Coma, which 
yield (Xcr) < 0.017 (Arlen et al. 2012), provided the CR uni- 
versality assumption holds. The most suitable cluster target 
for CR-induced ]/-ray emission, the Perseus cluster, cannot be 
used to competitively constrain the CR pressure using Fermi- 
LAT data since the Fermi-LP^ and MAGIC collaborations de- 
tected the central radio galaxy NGC 1275 in }/ rays in the energy 
range from 300 MeV to >300 GeV {Fermi: 300 MeV-300 GeV, 
MAGIC: >200 GeV) (Abdo et al. 2009; Aleksic et al. 2012b). 

Assuming both universality as well as the scaling relation we adopt 
throughout our work, which is characterized through Ay for the combined 
sample, for Coma, we find (Acr) < 0.011. 


Non-observations of y rays from the Perseus cluster above these 
energies by the MAGIC Collaboration (Aleksic et al. 2012a, 
2010b) provide limits similar to those obtained from analyzing 
LAT observations of Coma, {Xqv) < 0.017 alone (Arlen et al. 
2012; Zandanel & Ando 2013). 

Our limits on Xqr probe the entire ICM and are much more 
constraining in comparison to the limits on the non-thermal 
pressure contribution of the central ICM in several nearby cD 
galaxies that have been derived by comparing the gravitational 
potentials inferred from stellar and globular cluster kinematics 
and from assuming hydrostatic equilibrium of the X-ray emitting 
gas (Churazov et al. 2008, 2010). Depending on the adopted 
estimator of the optical velocity dispersion, they And a non- 
thermal pressure bias of Xnt = Pnt/Pth ~ 0.21-0.29, which 
probes the cumulative non-thermal pressure contributed by 
CRs, magnetic flelds, and unvirialized motions. It is, however, 
conceivable that those central regions of CC clusters, which 
probe the enrichment of non-thermal components as a result of 
AGN feedback (e.g., Pfrommer 2013), are characterized by a 
larger CR pressure contribution in comparison to the bulk of the 
ICM that probes CRs accelerated by shocks associated with the 
growth of structure and magnetic flelds that are reprocessed by 
these shocks. 

Complementary limits on CRs are also derived from ra- 
dio (synchrotron) observations (Pfrommer & EnBlin 2004a; 
Brunetti et al. 2007). Assuming a central cluster magnetic 
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Table 7 

Individual Flux Upper Limits (Alternative CR Profiles) 


Cluster 
Spatial Model 

a\JL 

ICM 

pCR 

y,exp 

ICM 

uUL 

^y,500MeV 

ICM 

TS 

ICM 

aCL 

Flat 

pCR 

y,exp 

Flat 

UUL 

^y,500MeV 

Flat 

TS 

Flat 

2A0335 

0.85 

1.8 

1.5 

0.0 

5.14 

0.4 

2.2 

0.0 

A0085 

1.64 

2.4 

3.9 

0.6 

6.15 

0.8 

4.6 

1.1 

A0119 

3.91 

1.3 

5.2 

2.7 

9.55 

0.6 

6.1 

3.0 

A0133 

3.83 

0.5 

2.0 

0.0 

15 

0.1 

1.9 

0.0 

A0262 

2.39 

1.0 

2.3 

0.0 

5.63 

0.7 

3.8 

0.0 

A0400 

59.23 

0.4 

22.5 

52.9 

95.1 

0.3 

24.5 

61.4 

A0478 

1.57 

1.9 

3.0 

0.0 

6.9 

0.4 

2.6 

0.0 

A0496 

2.74 

2.1 

5.7 

2.3 

14.66 

0.6 

9.0 

5.4 

A0548e 

12.8 

0.2 

2.7 

0.1 

17.91 

0.2 

2.7 

0.1 

A0576 

5.25 

0.6 

2.9 

0.3 

15.01 

0.2 

3.1 

0.1 

A0754 

3.8 

1.5 

5.6 

2.3 

19.61 

0.4 

7.1 

2.8 

A1060 

2.63 

1.7 

4.5 

0.4 

13.74 

0.6 

7.7 

0.1 

A1367 

3.68 

1.8 

6.4 

5.6 

18.64 

0.4 

7.7 

5.0 

A1644 

2.77 

1.2 

3.3 

0.0 

7.67 

0.5 

4.1 

0.0 

A1795 

0.58 

2.3 

1.3 

0.0 

5.03 

0.3 

1.6 

0.0 

A2065 

4.99 

0.8 

4.2 

1.5 

26.23 

0.2 

4.5 

1.5 

A2142 

0.57 

3.0 

1.7 

0.0 

2.61 

0.8 

2.1 

0.0 

A2199 

1.84 

2.4 

4.4 

1.9 

7.09 

0.7 

4.9 

1.2 

A2244 

2.23 

0.6 

1.4 

0.0 

6.86 

0.2 

1.4 

0.0 

A2255 

5.67 

0.8 

4.5 

5.1 

11.27 

0.4 

4.5 

4.4 

A2256 

0.49 

2.4 

1.2 

0.0 

2.13 

0.7 

1.6 

0.0 

A2589 

5.06 

0.5 

2.7 

0.0 

21.48 

0.2 

3.8 

0.0 

A2597 

1.97 

0.5 

1.0 

0.0 

12.69 

0.1 

1.1 

0.0 

A2634 

6.84 

0.6 

3.8 

0.4 

15.34 

0.3 

4.2 

0.3 

A2657 

3.01 

0.7 

2.1 

0.0 

27.51 

0.2 

4.1 

0.1 

A2734 

3.18 

0.4 

1.2 

0.0 

7.4 

0.2 

1.3 

0.0 

A2877 

2.14 

0.4 

0.9 

0.0 

16.41 

0.1 

1.2 

0.0 

A3112 

7.57 

0.8 

6.0 

13.8 

29.98 

0.2 

6.1 

13.4 

A3158 

1.51 

1.1 

1.6 

0.0 

3.78 

0.5 

1.7 

0.0 

A3266 

1.34 

2.9 

3.9 

2.1 

6.17 

0.7 

4.6 

1.7 

A3376 

8.43 

0.6 

5.4 

0.5 

28.36 

0.3 

7.5 

1.2 

A3822 

5.27 

0.4 

2.1 

0.0 

9.07 

0.2 

2.2 

0.0 

A3827 

0.72 

0.9 

0.6 

0.0 

3.91 

0.2 

0.7 

0.0 

A3921 

5.96 

0.4 

2.5 

0.6 

19.58 

0.1 

2.7 

0.6 

A4038 

2.49 

1.0 

2.4 

0.1 

11.23 

0.3 

3.6 

0.5 

A4059 

2.66 

0.7 

2.0 

0.1 

11.15 

0.2 

2.1 

0.0 

COMA 

0.41 

10.4 

4.3 

0.9 

1.37 

4.9 

6.8 

2.2 

EXO0422 

0.92 

0.6 

0.5 

0.0 

4.83 

0.1 

0.6 

0.0 

FORNAX 

7.19 

0.7 

4.7 

0.5 

74.75 

0.2 

12.3 

1.8 

HCG94 

8.44 

0.3 

2.6 

0.0 

23.26 

0.2 

3.7 

0.0 

HYDRA-A 

3.68 

1.2 

4.2 

3.0 

18.86 

0.2 

3.9 

1.2 

IIIZw54 

12.94 

0.4 

4.9 

0.4 

52.7 

0.1 

4.2 

0.0 

IIZwl08 

3.66 

0.4 

1.3 

0.0 

6.43 

0.2 

1.4 

0.0 

NGC 1550 

5.38 

0.5 

2.5 

0.0 

25.85 

0.2 

4.9 

0.0 

NGC 5044 

9.95 

0.3 

3.1 

0.0 

61.79 

0.1 

6.1 

0.1 

RXJ2344 

5.79 

0.5 

2.9 

1.3 

30.83 

0.1 

2.8 

0.4 

S405 

3.5 

0.4 

1.3 

0.0 

6.16 

0.2 

1.5 

0.0 

S540 

14.76 

0.3 

4.1 

1.2 

48.07 

0.1 

4.5 

1.0 

UGC03957 

10.74 

0.4 

4.1 

0.5 

75.34 

0.1 

5.3 

0.6 

ZwC11742 

2.34 

1.0 

2.3 

0.0 

12.66 

0.2 

3.0 

0.0 


Notes. Individual flux predictions and upper limits for alternative CR profiles 
for photon energies above 500 MeV. Columns 2-5 refer to the individual scale 
factor, the flux prediction above 500 MeV (F^exp)’ derived upper limit 
on the y-ray flux, and the individual TS value for the ICM model, while 
Columns 6-9 represent the values obtained for the flat CR model. All photon 
fluxes are given in 10“^^ photon s“^ cm“^. 


field of 5 ~ 1 /xG and CR spectral indices a = 2. 1-2.4, 
this approach allows the CR energy to be constrained to a few 
percent while the limits are less stringent for steeper spectra and 
lower magnetic fields (Aharonian et al. 2009a). 


5.(5. Hadronic Injection Efficiency 

The distributions of CRs within the virial regions of clusters 
are built up from shocks during the cluster assembly (Ryu 
et al. 2003; Pfrommer et al. 2006; Pinzke et al. 2013). While 
strong shocks are responsible for the high-energy population 
of CRs that could potentially be visible at TeV y-my energies, 
intermediate Mach-number shocks with M — 3-4 build up the 
CR population at GeV energies (Pinzke & Pfrommer 2010), 
which can be constrained through Fermi observations. The 
normalization of both the CR population and the ]/-ray flux 
scale with the acceleration efficiency at shocks of corresponding 
strength. Our fiducial CR model (Pinzke & Pfrommer 2010) uses 
a simplified model for CR acceleration (EnBlin et al. 2007) in 
which the efficiency rises steeply with Mach number for weak 
shocks and saturates already at shock Mach numbers > 3. 
Observations of supernova remnants (Helder et al. 2009) and 
theoretical studies (Jones & Kang 2005) suggest a value for 
the saturated acceleration efficiency of ^p,max — 0-5, i.e., the 
fraction of shock-dissipated energy that is deposited in CR ions. 
Following these optimistic predictions, the model of Pinzke & 
Pfrommer (2010) assumes this value for the saturated efficiency, 
which serves as an input to our analysis. This model provides 
a plausible upper limit for the CR contribution from structure 
formation shocks in galaxy clusters since more elaborate models 
of CR acceleration predict even lower efficiencies for the Mach- 
number range M. — 3-4 of relevance here (Kang & Ryu 2013, 
see also Appendix C for a more detailed discussion). If we 
scale the CR pressure contribution linearly with the maximum 
acceleration efficiency, using the previously derived limit on 
Ay, we find combined sample while 

the CC- and NCC-subsamples yield 25% and 24% maximum 
acceleration efficiency, respectively. We note that this constrains 
the maximum acceleration efficiency only in the simplified 
model adopted here. For different acceleration models, these 
upper limits provide conservative constraints on the acceleration 
efficiency at intermediate strength shocks of Mach number 
M — 3-4, a regime complementary to that studied at supernova 
remnant shocks. 

However, these conclusions rely on two major assumptions, 
namely, CR universality and the absence of efficient CR trans- 
port relative to the plasma rest frame. The latter assumption 
hypothesizes that CRs are tied to the gas via small-scale tangled 
magnetic fields, which implies that they are only advectively 
transported and that we can neglect the CR streaming and dif- 
fusive transport relative to the rest frame of the gas. Early work 
on this topic suggests that such CR transport processes are at 
work in clusters and cause a flattening of radial CR profiles that 
can significantly reduce the radio and }/-ray emission at high 
energies probed by Cherenkov telescopes but remains largely 
unaffected at lower energies probed by Fermi (Wiener et al. 
2013). Moreover, different formation histories of clusters cause 
spatial variations of the CR distribution and hence a deviation 
from a universal distribution (Pinzke & Pfrommer 2010). To 
date, there is no consensus about the size of these effects, and 
more work is needed to fully quantify them. 

6. SYSTEMATIC UNCERTAINTIES 
6.1. Choice of Fiducials: Binning, Region Size, Free Sources 
6.1.1. Binning 

While the number of spectral bins is the same for the whole 
sample (nominally 18 bins), the number of spatial bins varies 
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Figure 11. Individual 95% upper limits on Xcr for each of our 50 galaxy clusters in this analysis assuming the jointly derived scale factor we obtain in our analysis 
for the full sample (blue diamond), CC clusters (red downward triangle), and NCC clusters (green circle). The dashed lines represent the median upper limit for the 
combined (blue) while the median upper limits for CC and NCC are the same (shown in black). 

(A color version of this figure is available in the online journal.) 


with the ROI size, because the bin width is constant (0?1). 
In addition, the extended templates used to model the cluster 
emission are also binned. To check the effect of binning, we vary 
the nominal values by 50%. Aside from potentially increasing 
the computation time for larger numbers of spatial bins, we 
find that our choice of binning does not change the results by 
more than ~1%. Similarly, we find that varying the number of 
spectral bins changes the resulting limits on A by at most 5%. 

6.1.2. Region Size 

We use ROIs of varying sizes, ranging from 8°-16° in radius 
(Table 2). To make sure that this choice does not introduce 
any significant bias, we compare the fitted values for all free 
parameters for the null hypothesis in these smaller regions 
with ROIs with 25% larger radii and find variations of these 
values which are less than 3% with respect to larger regions. 
However, we note that larger regions allow a more stringent 
determination of the background model which is reflected 
by smaller uncertainties on the Galactic and isotropic diffuse 
components than for the case of smaller regions (compare to 
error bars in Figure 5). 

6.1.3. Free Sources 

We choose to use the two year source list to model the 
data collected during four years of LAT observations. While 
we ensure that residual excesses are mitigated by allowing the 


normalizations of the known point sources to vary, the choice of 
leaving the normalizations of sources within 4° of each cluster to 
vary freely is somewhat arbitrary. Freeing the normalizations of 
only those sources within ^200 + 1° does not change our results 
on A by more than 10%. 

6.2. Event Classes and Instrument Response Functions 

The IRFs consist of three separate parts (see Ackermann et al. 
2012a, for details): the effective area which has an associated 
uncertainty of at most ~10% in the energy range we consider, 
the PSF whose uncertainty can be conservatively estimated 
to be ~ 15%, and the energy dispersion (whose uncertainty is 
negligible for this analysis). Using bracketing IRFs (see Sections 
5.7.1 and 6.5.1 in Ackermann et al. 2012a) to quantify these 
uncertainties, we find that while individual ROIs may show 
variations of up to ~21%, the effect on the combined scale 
factors and quantities derived from it is less than 7%. 

6.3. Dijfuse Emission 

Spatial residuals due to mismodeling of the large-scale 
Galactic diffuse foreground emission may be misinterpreted 
in terms of an extended ]/-ray excess. We compare results 
derived using the standard diffuse emission model adopted for 
the baseline analysis (based on empirical fits of multiple spatial 
templates to y-ray data) to results obtained when using a set of 
eight alternative diffuse emission models that were created using 
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Table 8 

Alternative Models for Galactic Diffuse Emission 


Label 

CR Source Distribution 

Halo Size 
(kpc) 

Spin Temperature 
(K) 

A 

Lorimer 

10 

105 

B 

Lorimer 

10 

150 

C 

Lorimer 

4 

1Q5 

D 

Lorimer 

4 

150 

E 

SNR 

10 

105 

F 

SNR 

10 

150 

G 

SNR 

4 

105 

H 

SNR 

4 

150 


Notes. Overview of the alternative diffuse models used for assessing the 
systematic uncertainties in the model for the Galactic diffuse emission. 
We chose to vary the three most important input parameters that were 
found in scanning the parameter space in Ackermann et al. (2012c). 


a different methodology with respect to the standard diffuse 
emission model 

We chose these models to represent the most important 
parameters scanned in Ackermann et al. (2012c), in particular, 
CR source distribution, halo size, and spin temperature. We 
summarize the properties of the alternative models in Table 8, 
and we refer readers to de Palma et al. (2013) for details. The 
models we employed were tuned to the P7REP data. Although 
the models were created such that different components along 
the 1.0. s. could be fit separately, we only adopted a free overall 
normalization, since at high Galactic latitudes the vast majority 
of the gas resides in the neighborhood of the solar system. 
Moreover, having different components as additional degrees of 
freedom in the fit makes a comparison of the TS values with the 
baseline analysis more difficult. 

We emphasize that these eight models do not span the com- 
plete uncertainty of the systematics involved with interstellar 
emission modeling. They do not even encompass the full un- 
certainty in the input parameters that are varied. The resulting 
uncertainty should therefore only be considered as one indicator 
of the systematic uncertainty due to interstellar emission mod- 
eling. The tests we performed using these additional models 
considered a different energy range, from 500 MeV-100 GeV, 
because the alternative models were not derived for higher en- 
ergies. However, since events at low energies dominate the fit, 
this difference is negligible. We have explicitly verified this by 
repeating the baseline analysis up to 100 GeV and we found that 
the differences between the computed combined scale factors 
are <1%. 

In Figure 12, we show how the combined upper limit on the 
scale factor, varies for the alternative models with respect 
to our standard model as well as how the significance of the 
excesses observed in A1367 and A3 112, and A400 changes 
when using the alternative models. 

The spread in limits for the scale factor is rather small 
between the alternative diffuse models, but the choice of using 
the standard diffuse model versus any of the alternative diffuse 
models can affect the inferred limits for the different cluster 
samples by 20%-30%. Comparing the TS values for the three 
clusters indicates small variations across the alternative models 
for A400, A3 112, and A1367. We repeated this procedure for 
the derivation of the individual flux limits (see Section 5.4) and 
found that the majority of our clusters follow the same trend. We 


We also use a different set of isotropic diffuse templates that were created 
in conjunction with the alternative Galactic diffuse templates. 
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Figure 12. In the upper panel, we show the 95% combined upper limit on Ay 
for the respective (sub)samples for the alternative diffuse models (A-H, see 
Table 8 and the accompanying text for details). In the bottom panel, we show 
the associated TS values for the three clusters that exhibit significant excess 
emission. 
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have marked the clusters which show variations beyond 30% in 
Table 6. 

We summarize the systematic uncertainties discussed in 
this section in Table 9, and we note that the main source of 
uncertainty is the accurate modeling of the foreground Galactic 
diffuse emission. 

7. SUMMARY AND CONCLUSIONS 

We have used a data set of 4yr of all- sky data from the 
Fermi-LAI detector and performed a search for high-energy 
]/-ray emission originating from 50 X-ray luminous galaxy 
clusters. We specifically consider hadronically induced y rays 
originating from the ICM as described by the universal CR 
model by Pinzke & Pfrommer (2010) and employ a joint 
likelihood analysis to constrain the normalization of a common 
scale factor among clusters that is theoretically expected to 
describe the y-my luminosity of the ICM. In order to allow 
for different emission scenarios, we categorize clusters in our 
sample by their morphologies and separately consider CC and 
NCC subsamples. 

We find evidence for excess emission at a significance of 
2.7a, naively taken as a first indication of y-ray emission from 
galaxy clusters. However, upon closer investigation, we find 
that this global significance originates mainly from individual 
excesses present in three galaxy clusters: A1367, A3 112, and 
A400, the latter yielding a post-trial significance of 6.7a alone. 
For these three clusters, the LAT data alone cannot conclusively 
support or reject the hypothesis that the excess emission arises 
from within the clusters. The best- fit location of each excess is 
located within the virial radius of the respective cluster, but it is 
offset from the cluster center. 

With respect to the universal CR model, we also note that 
the associated scale factors are significantly larger than the ones 
derived from other clusters in the sample. We also argue that in 
all three clusters there are individual radio galaxies which may 
be the origin of observed excesses. 

We establish bounds on the common scale factor Ay, and 
we use this to derive individual upper limits on the y-ray flux. 
In addition, we use the jointly derived limit on to calculate 
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Table 9 

Systematic Uncertainties 
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Type 

Variation of Input Parameters 

Impact on Results 

Spectral bins 

±50% 

<5% 

Spatial bins 

±50% 

<1% 

Spatial template bins 

±50% 

<1% 

Small ROIs 

+25% 

~3% 

Number of free sources 

4° ^ 0200 + 

<10% 

IRF uncertainties: effective area 

±10%^ 

<7% 

IRF uncertainties: PSF 

±15%’^ 

<4% 

Diffuse model uncertainties 

Alternative diffuse models^ 

15%-25% more stringent limits 


Notes. Overview of systematic uncertainties as discussed in Section 6. We note that the largest impact on 
the results is due to the model for the Galactic diffuse emission. 

^ We chose a radius of 4° around each cluster center to account for photon contamination due to the PSF 
at low energies. In this test, we modify the radius in which we leave the normalization free to vary within 

O200 + 1 °- 

^ We employ the bracketing IRF approach as discussed in Ackermann et al. (2012a) and use the tabulated 
values to scale the relevant IRF components. 

^ We use a set of alternative diffuse emission models and replace the standard emission template used in 
the baseline analysis with these. 


limits on the volume- averaged CR-to-thermal pressure (Xcr). 
We compute median upper limits calculated within R 200 , with 
the most stringent one being (Xcr) < 0.012 for the combined 
sample and (Xcr) < 0.013 and (Xcr) < 0.014 for the CC and 
NCC sub-samples, respectively. Assuming a linear dependence, 
our limits on translate into a combined limit of the hadronic 
injection efficiency, Cp,inj. by large-scale structure formation 
shocks in the Mach number range M — 3-4 to be below 
21% for the combined sample and 25% and 24% for the CC 
and NCC clusters, respectively. Removing the aforementioned 
three clusters that exhibit excess emission provides even more 
stringent limits on Ay, which for the combined sample, yields 
AJ/^ = 0.29 that translates into (Xcr) < 0.008 within R 200 
and ^p,inj(A^ - 3-4) < 15%. Our limits on (Xcr) and ^p,inj 
are the most stringent to date, constraining hadronic emission 
scenarios that predict astrophysical y rays originating in the 
ICM of galaxy clusters. 

The systematic uncertainty associated with the modeling of 
the Galactic foreground emission represents the largest source 
of uncertainty affecting limits on extended emission from the 
ICM presented in this work. To account for this, we have tested 
our results against a set of alternative diffuse models spanning a 
range of interstellar emission model parameters. We find that the 
alternative models provide limits that differ from the baseline 
analysis by 20%-30%. 
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APPENDIX A 

MC SIMULATION STUDIES 

We use the simulation package gtobssim to efficiently gen- 
erate MC realizations of the y-ray sky using the parametrized 
instrument response. 

A. 7. Minimum Energy Threshold 

The joint likelihood approach discussed in Section 3.1 makes 
the assumption that the individual data samples are uncorrelated. 
Assuming two sources and ^’ 2 , in the uncorrelated case, the 
composite /7-value is /?x = P\ ^ Pi, where p\ and p 2 are 
the /7-values associated with s\ and 5 * 2 , respectively. This case 
corresponds to two sources that are far away from each another. 
In this case, the difference between /7joint, which is the derived 
/7-values from the joint likelihood, should be minimal, while for 
small distances, correlations impact the derivation of /7joint, and 
thus lead to an overestimation of the significance associated with 
this /7-value. We assess this bias using an isotropic simulation 
with two identical power-law sources with E = —2.3 that 
were each modeled as an extended source with a disk of 2° 
in diameter. In Pigure 13, we show the difference between 
/7joint and pj^. We find that at small distances (<3°), for all 
minimum energy thresholds, there is a substantial bias due to 
overlaps. Toward larger distances, this bias is reduced. In this 
toy model, 3° corresponds to 7?2oo + 1°- However, by requiring 
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Figure 13. We show the difference between the jointly derived p-value and 
the case of uncorrelated samples, where pi: = p\ pi (see text for details) 
simulating two power-law sources with identical spectral model and modeled as 
a disk with 2° in diameter for a minimum energy threshold of (100 MeV, 
200 MeV, 500 MeV, and 1 GeV). The solid line corresponds to our analysis 
threshold, chosen to minimize this overlap bias while maximizing the expected 
y-ray flux. 


larger distances between sources, we reduce the number of 
viable cluster candidates for our search. Hence we decided to 
use £min = 500 MeV as this threshold maximizes the number 
of clusters to be included (and thus the expected signal) while 
minimizing the bias on j^joint- 

A.2. Significance Assessment 

In Figure 14, we show ROI-specific TS-distributions for 
a background only simulation (top). We find that the 
TS -distribution in the background-only case can be well de- 
scribed by (1/2)8 -f (l/2)Xk and obtain k = 1.1 zb 0.1 in an 
unbinned maximum likelihood fit, giving rise to the usual defi- 
nition of significance as VtS. 

A. 3. Signal Studies 

In addition, we include the results from simulations including 
a (weak) putative CR-induced y-ray signal corresponding to 
Ay = 0.29 (nominal best-fit value for the combined sample) 
as well SiS Ay = 1 (middle panel in Figure 14), i.e., assuming 
the predictions by Pinzke et al. (2011). Finally, we repeat the 
analysis assuming a strong signal, characterized by = 10 
(bottom panel in Figure 14). For all simulations, the true value of 
Ay is recovered in the combined result, validating our analysis 



TS TS 


Figure 14. In the top left panel, we show the TS-distribution for flve background- only simulations and fit the distribution to the null hypothesis according to our 
analysis (refer to Section 3.1 for details). In the top right panel, we show the constraints on the y-ray flux using the ROI-specific individual scale factors Ayj and 
relate them to the associated TS values, denotes the Spearman rank correlation coefficient and sig refers to the combined significance derived from TSgiobai- The 
middle and lower panel show the same but assume a putative y-ray source in addition to the background. We show that for = 1 and even more so for A ^ = 10 the 
TS-distributions clearly depart from a distribution. It should be noted, however, that while the background simulation is based on flve MC realizations of the y-ray 
sky, the various signal simulations represent a single example realization for each assumed signal. 
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approach, although the global significance varies with respect to 
the signal simulation. For the background-only case, only upper 
limits can be derived. Assuming the nominal best-fit value for 
Ay from the combined sample, yields a combined significance 
of2.1a. For = landA^^ = 10, these values are much higher, 
yielding a global TS value corresponding to a significance of 
13.3cr and 19.9a, respectively. 

Given that a simulation with Ay = 0.29 yields a combined 
significance comparable with what we have found in our analy- 
sis, we investigated how this signal could be studied further. To 
that end, we compare the expected y-ray flux based on the ROI- 
specific scale factors Ay , i with their associated ROI TS val- 
ues and quantify the correlation through calculating Spearman- 
rank correlation coefficient (Spearman 1904), denoted by r^. 
We find Spearman-rank coefficients >0.5 indicating a corre- 
lation. However, even for the background case, = 0.6 is 
obtained, which is the same as what we find in the signal 
case for Aj, = 1. This illustrates that a correlation analysis 
with such a weak signal, or further, studying the excess we 
find, is difficult. Only for the strong signal of A^^ = 10 do 
we find a correlation coefficient = 0.94 indicating a strong 
correlation between expected ]/-ray flux and associated ROI 
TS value. 


APPENDIX B 

ANALYTIC COSMIC-RAY MODEL 

Eollowing the analytic CR formalism (Pfrommer & EnBlin 
2004a; Pfrommer et al. 2008; Pinzke & Pfrommer 2010), we 
obtain the volume- weighted, energy-integrated, and omnidirec- 
tional (i.e., integrated over the 47T solid angle) y-ray source 
function due to pion decay, 

Xj^Q_y{R, E) = Kj^Q_y{R)Xj^0_y{>E). (B 1 ) 

Here the spatial part of the y-ray emission is determined by 

Kj^Q_y{R) = CcR(R)Pgas(^), (B2) 

where the CR proton distribution is 


Pgas(^) y V^trans/ J 

-1 

C'center ? and 

(B3) 

C200 = 1-7 X 10-’ X (M 2 oo/ 10 '^ 

(B4) 

^trans = 0.021 /?200 X (M20o/10‘^ 

(B5) 

= 1.04 X (M2oo/10‘^ 

(B6) 


where C = Cnipl p denotes the dimensionless normalization 
of the CR distribution function. Eor massive clusters (M 200 ~ 
10^^ Mq) the CR distribution traces the gas density, while the 
CR density is slightly enhanced in the center for smaller systems. 
Note, however, that the y-ray flux depends only weakly on the 
exact CR density in the center (i.e., Ccenter) since most of the 
flux originates from outside the transition region. 


The spectral part of Equation (Bl) for the photon energies 
relevant to Eermi-hAI (lOOMeV ^ Ey < I TeV) is given by 


Xj,0_y(>E) 



(B7) 


where the sum over i extends over A = (0.767, 0.143, 0.0975), 
and a = (2.55,2.3,2.15). The y-ray spectrum rises until a 
maximum at the mass of the meson followed by a concave 
shaped tail determined by the universal CR spectrum. The shape 
parameter 8t ^ 0. 14 + 0.44 allows us to accurately predict 

the emission close to the pion bump in combination with the 
effective inelastic cross-section for proton-proton interactions, 
cTpp / 32(0.96 - 1 - mbam. We have also introduced 

the equation 


[B, (a, bXl = B,, {a, b) - B,, (a, b) , (B8) 

where Bx {a, b) denotes the incomplete Beta- function. 

APPENDIX C 

RELATIONSHIP BETWEEN CR ACCELERATION 
EEEICIENCY AND GAMMA-RAY ELUX 

In this Appendix, we investigate how we can use upper 
limits on the y-ray flux to constrain the CR injection effi- 
ciency in various models for CR acceleration. The CR model 
adopted in this paper (Pinzke & Pfrommer 2010) is based 
on a simplified scheme (EnBlin et al. 2007) to compute the 
CR-energy acceleration efficiency at shocks (in units of the 
shock-dissipated thermal energy, corrected for adiabatic com- 
pression), ^p,inj(7W) = £cR/^diss- H cmploys the thermal leak- 
age model (e.g., Ellison & Eichler 1984; Berezhko et al. 1994; 
Kang & Jones 1995), which conjectures a momentum thresh- 
old for injection that is a constant multiple (xinj = 3.5) of the 
peak thermal momentum, at which the CR power-law distri- 
bution connects to the post- shock Maxwellian. More refined 
models (such as in Kang & Ryu 2011, 2013) are motivated 
by nonlinear shock acceleration, and they fix the injection mo- 
mentum by the considering that the particle speed should be 
several times larger than the downstream flow speed in or- 
der for suprathermal particles to diffuse upstream across the 
shock transition layer. This yields a Mach-number dependent 
Xinj, which increases for weaker shocks such that a progres- 
sively smaller fraction of particles can participate in the process 
of diffusive shock acceleration. As a result, for the preferred 
values of the ratio of the downstream turbulent-to-background 
field, sb ^ 0.25, those models predict ^p,max ^ 0.4-0. 5, de- 
pending on the existence of a pre-existing CR population.^^ 
While those values for the acceleration efficiency are now chal- 
lenged by y-ray observations of supernova remnants, additional 
physics (such as amplification mechanisms of the magnetic 
field) may lower the value oi Sb and cause the acceleration 
efficiency to saturate at lower values (Kang & Ryu 2013). 

This is obtained by taking the ratio of CR acceleration efficiency {rf) and 
gas thermalization efficiency (S) in the limit of large Mach numbers (Kang & 
Ryu 2013). 
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^inj 



Figure 15. Acceleration efficiency ^p,mj as a function of shock Mach number 
(Ai) for two different post- shock temperatures of 0.1 keV (thick bright) and 
1 keV (thin faded). We show the acceleration efficiency that was used in our 
simulations (^p,max = 0.5, red solid) and that was scaled to different values for 
the saturated acceleration efficiency (^p,max = 0.1 and ^p,max = 0.05; red dotted 
and dashed, respectively). Also shown in blue and green is the acceleration 
efficiency in a model by Kang & Ryu (2013) for different values of the ratio 
of downstream turbulent-to-background magnetic held, which determines the 
injection efficiency in their model. Since shocks with Mach numbers A4 — 3-4 
are mostly responsible for injecting CRs in clusters, constraints on ^p^nj in the 
model by EnBlin et al. (2007) are more conservative in comparison to the model 
by Kang & Ryu (2013). 

(A color version of this hgure is available in the online journal.) 


Moreover, in weak heliospheric shocks, additional shock phe- 
nomena (e.g., whistler waves in the shock front, etc.) are ob- 
served that may add to the uncertainty of the acceleration ef- 
ficiency. In summary, while the models for the acceleration 
efficiency in shocks became more detailed and physical over 
the last year, new observations point to the necessity of further 
improving those models. 


In Figure 15, we show the relation between acceleration 
efficiency and shock Mach number. For the relevant energy 
regime that we consider in this work (Ey > 500 Me V), our CR 
model predicts a CR spectral index of ~2.4 that flattens toward 
higher energies (Pinzke & Pfrommer 2010). This implies that 
shocks with Mach numbers Ai ^3.5 (depending somewhat 
on the post-shock temperature) are the most relevant for the 
CR budget. Because the acceleration efficiency has already 
saturated for this range of shock strengths in the model of 
EnBlin et al. (2007), the CR pressure and thus the hadronically 
induced y-ray luminosity are approximately proportional to 
^p,max- For different acceleration models (such as in Kang & 
Ryu 2013), these upper limits provide interesting constraints 
on the acceleration efficiency at intermediate strength shocks 
of Mach number A4 ~ 3-4, a regime complementary to that 
studied at supernova remnant shocks. 

APPENDIX D 

RADIO AND X-RAY SOURCES IN THE FIELD 
OE VIEW OE A3 112, A1367 AND A400 

In this section, we provide a discussion of the radio sources in 
the field of view of the three clusters as discussed in Section 5.3. 
We note that a detailed characterization of the origin of the 
excess emission is beyond the scope of this work. Based on the 
refined best-fit positions from our higher resolution TS-map, we 
performed a search for sources within the 3cr contours as shown 
in Eigure 16. Below, we provide supplemental information 
regarding the three clusters discussed in the main text. 

1. Eor A 1367, our best-fit position is consistent with that of 
the radio galaxy 3C264 (Eey et al. 2004). As we cannot 
distinguish between A1367 and 3C264, we conclude by 
similarity with previous y-ray detections in clusters that 
we likely observe y-ray emission from the radio galaxy 
(e.g., M87 in Virgo or IC 310 in Perseus). 

2. Similarly, the origin of the emission toward A3 112 may 
be from the radio galaxy PKS 0316-444, which is located 
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Figure 16. TS maps of our three cluster candidates with TS > 9. Each map shows a 0?5 x 0?5 section (0?6 x 0?6 for A400) that was recentered to the best- fit position 
obtained using gtf indsrc. The best-fit position from the refined TS calculations is marked as red x. Shown in red are the 3cr (solid) and 4cr contours (dashed). The 
blue diamond-shaped points correspond to the NED positions of the radio galaxies as discussed in the text. The purple upright open triangles denote Chandra X-ray 
sources that fall within the error circle of the y-ray point source. Overlaid radio contours (blue dotted) were obtained from the NVSS for A1367 and A400 (Condon 
et al. 1998). The radio contours for A3 112 were obtained from the Parkes-MIT-NRAO Survey (Condon et al. 1994). The cluster center is marked by a black cross in 
each panel. The virial radius of the cluster is indicated by the dashed black line (partially visible in the maps for both A31 12 and A400). Coordinates are taken from 
NED/SIMBAD. Each pixel is 0?2 across. 

(A color version of this figure is available in the online journal.) 
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Figure 17. Fitted counts spectra for 9 of the 26 analysis ROIs without additional cluster sources. The crosses indicate the measured counts in each energy bin while 
the black dashed lines show the total sum of all model counts for all components. The gray dashed lines refer to the Galactic diffuse component while the gray dotted 
lines correspond to the isotropic extragalactic diffuse component. Solid lines indicate additional background sources. We obtain reasonable fits in all energy bins. The 
lower panel for each ROI shows fractional residual counts integrated over the entire ROI. 


4'.2 away from the best-fit position of the excess (Costa & 
Loyola 1998). 

3. In the vicinity of the best-fit position of the excess to- 
ward A400, a (not further classified) radio source NVSS 
J025857-f055240 was reported by Condon et al. (1998). 
Because of this positional coincidence, it is plausible to at- 
tribute the observed ]/-ray emission to this object, although 
this hypothesis warrants further investigation. 

In addition to the previously discussed radio sources, there are 
multiple Chandra X-ray sources that fall within the error circle 
of the best-fit positions for the excesses in A1367 and A3112, 
respectively. While the association of the excesses in y rays 


with individual radio galaxies is well in line with previous 
y-ray detections, e.g., M87 in Virgo or NGC 1275 and IC 310 
in Perseus, the similarity argument we present here is not 
sufficient to claim detection of y rays from these respective 
objects. 

APPENDIX E 

ROI-SPECIEIC COUNTS/MODEL COMPARISON 

We provide for each ROI the observed photon counts 
in each energy bin along with the predicted model counts 
from the best-fit background-only model and show this in 
Eigures 17-19. 
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Figure 18. Fitted spectra for ROIs 10-18. See caption of Figure 17 for details. 
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Figure 19. Fitted spectra for ROIs 19-26. See caption of Figure 17 for details. 
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